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Abstract 


A new frequency-domain system identification algorithm is presented 
for system identification of substructures, such as payloads to be flown aboard 
the Space Shuttle. In the vibration test, all interface degrees of freedom where 
the substructure is connected to the carrier structure are either subjected to 
active excitation or are supported by a test stand with the reaction forces 
measured. The measured frequency-response data is used to obtain a linear, 
viscous-damped model with all interface-degree of freedom entries included. 
This model can then be used to validate analytical substructure models. This 
procedure makes it possible to obtain not only the fixed-interface modal data 
associated with a Craig-Bampton substructure model, but also the data associ- 
ated with constraint modes. With this proposed algorithm, multiple-boundary- 
condition tests are not required, and test-stand dynamics is accounted for 
without requiring a separate modal test or finite element modeling of the test 
stand. Numerical simulations are used in examining the algorithm’s ability 
to estimate valid reduced-order structural models. The algorithm’s perfor- 
mance when frequency-response data covering narrow and broad frequency 
bandwidths is used as input is explored. It’s performance when noise is added 
to the frequency-response data and the use of different least squares solution 
techniques are also examined. The identified reduced-order models are also 
compared for accuracy with other test-analysis models and a formulation for a 
Craig-Bampton test-analysis model is also presented. 
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Chapter 1 


Introduction 


The dynamic behavior of complex engineering structures can be pre- 
dicted by using various analytical and numerical methods, with the most pop- 
ular method being the finite element method. A finite element model (FEM) 
is an approximation of a continuous structure, and assumptions are made in 
the development of the FEM concerning the distribution of mass and stiffness, 
approximation of boundary conditions, estimation of material properties, and 
so forth. As a result, the FEM may or may not accurately represent the original 
structure. Before this mathematical model can be used to perform any sort of 
analysis and the results used with confidence, the model must be test- validated 
to ensure that it accurately represents the physical structure. Thus, there ex- 
ists the need to be able to characterize the dynamic behavior of structures 
experimentally. 

For some structures, performing a vibration test on the entire struc- 
ture is usually not feasible or practical due to the size and complexity of the 
structure. To facilitate testing, the original structure is divided into smaller, 
more manageable components, referred to as substructures. Component mode 
synthesis, or substructuring, techniques are generally employed to analyze the 
individual components and couple them together to form an analytical model 
of the original structure [1, 2]. For example, the Space Shuttle Orbiter, pay- 
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load, external tank, and solid rocket boosters are each substructures which can 
be analyzed individually and coupled together to form the final launch vehicle 
model. 

To verify an individual substructure model, a modal test is performed 
to identify the modal parameters of the actual substructure. The modal pa- 
rameters of interest are the natural frequencies, damping factors, natural mode 
shapes, and the residual mass and stiffness. Since the natural frequencies 
are scalar quantities, they can be readily compared. The comparison of the 
mode shapes is not as straightforward and involves orthogonality and cross- 
orthogonality checks using the experimental mode shapes and test-analysis 
models. 

When the vibration test is performed, the boundary conditions should 
be representative of those experienced during actual flight conditions. It is im- 
portant that that the correct boundary conditions be provided at the substruc- 
ture interface locations in order to obtain accurate results. For Space Shuttle 
payloads, this has sometimes led to fixed-base modal tests being performed to 
verify the analytical substructure model. The difficulty with this test method 
is that it requires the interface locations to have zero translational and rota- 
tional displacement. In actual test conditions, this situation is very difficult 
to achieve, since the test stand used to support the test article is rarely rigid 
enough and the dynamics of the test stand must be considered. 

Therefore, to eliminate these problems, a new test method and al- 
gorithm called the Substructure System Identification (SSID) algorithm has 
been developed to identify substructures and verify analytical substructure 
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models [3, 4]. In the vibration test, all interface degrees of freedom are either 
actively excited by a shaker or connected to a test stand with the reaction forces 
measured. The results can be used to obtain a linear, visoous-damped model of 
the substructure. The purpose of this investigation was to verify the proposed 
algorithm’s ability to identify valid structural models. Numerical simulations 
were used to simulate test data, and the algorithm’s performance when noise 
was added to the test data, and the effects of spatial and frequency truncation, 
were examined. 

Chapter 2 reviews pertinent test methods and structural system iden- 
tification algorithms. In Chapter 3, the theory of the Substructure System 
Identification algorithm is discussed. Model reduction and various test-analysis 
models used to aid in the comparison between the test and analytical results 
are discussed in Chapter 4. Chapter 5 highlights some of the computational 
aspects of the SSID algorithm that affect the final solution, such as different 
methods for solving an over-determined system of linear equations. In Chap- 
ter 6, the results of various numerical simulations using the SSID algorithm 
are presented and analyzed. Finally, in Chapter 7, conclusions are drawn and 
recommendations for future research are proposed. 



Chapter 2 


Review of Literature 


Reviews of previous test methods and structural system identifica- 
tion algorithms that are relevant to this work are given in this chapter. Test 
verification methods are discussed in Section 2.1. Section 2.2 reviews several 
different time-domain and frequency-domain structural system identification 
algorithms. 

2.1 Test Verification of Analytical Substructure Mod- 
els 

2.1.1 Test Verification of Finite Element Structural Models 

The response of complex structures subjected to dynamic loading has 
been a subject of interest for many decades. The finite element method is 
usually employed to facilitate the analysis of the structure. A finite element 
model of the continuous structure is created by discretizing it to A finite degrees 
of freedom (DOF). The resulting second-order differential equations of motion 
for a linear time-invariant model can be written in the form 

Mx{t ) + Cx(t) + Kx(t) = Df(t) (2.1) 

where M,C, and K are respectively the mass, damping, and stiffness matrices 
of the system, x(t) is the time-dependent coordinate vector, D is the force- 
distribution matrix, and f(t) is the time-varying applied force. Based on the 
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following physical assumptions about the structure, the system matrices are 
assumed to be positive definite or semidefinite: 1) all degrees of freedom have 
inertia, so the system kinetic energy is strictly positive and the mass matrix is 
positive definite, 2) for any forced motion, the system does not create energy, so 
the damping matrix is positive semidefinite, and 3) the strain energy is always 
positive, so the stiffness matrix is positive semidefinite. 

The FEM must provide a physically significant representation of the 
physical structure. However, assumptions are made in the development of the 
FEM, and much expertise is required in order to obtain reliable and valid 
results. The initial FEM must be experimentally validated before it can be 
used to predict the response of the structure. 

An experimental model can be obtained using one of the identification 
algorithms described in the next section. Then, a comparison between the 
experimental model and an analytical test-analysis model (TAM) is made in an 
attempt to verify the math model. Test-analysis models are discussed further 
in Chapter 4. 

The analytical model is deemed to be a good representation of the 
physical structure if there is good agreement between the measured and an- 
alytically computed modal parameters. Generally, if the identified natural 
frequencies are within 10%, and an orthogonality check of the test modes with 
respect to a TAM mass matrix is satisfied, the math model is considered “test- 
verified.” Acceptable orthogonality results are characterized by values of all of 
the off-diagonal terms being less than 0.1 when the diagonal values are normal- 
ized to unity [5]. The test data is generally considered to be correct, therefore 
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any discrepancies between the analytical and the measured modal parameters 
are attributed to the math model. 

Numerous procedures exist for updating the finite element model to 
match the experimental data. There are two general approaches for model 
updating. One approach is to modify the entire matrices, and the other is 
to update only individual elements within the matrices [6, 7]. The analytical 
mass matrix is generally considered to be a good estimation of the structure s 
mass characteristics, so many updating techniques focus on updating only the 
stiffness matrix. There also exist numerous methodologies that update the 
analyti cal model so as to minimize the difference between the analytical and 
experimental frequency response functions [8]. 

2.1.2 Space Shuttle Payload Model Verification Tests 

A coupled-loads analysis for lift-off and landing events must be per- 
formed for any payload that will fly aboard the Space Shuttle. The mathe- 
matical model used to perform the coupled-loads analysis must be in the form 
of Craig-Bampton modes or mass and stiffness matrices [9]. Before using the 
math model in the coupled-loads analysis, the math model must be test- verified 
by correlating it with modal test data. 

Test-verification of math models has traditionally been accomplished 
by performing a modal test on the payload with the interface degree of freedom 
constrained. For a proper validation, the boundary conditions imposed during 
the test should match those in service, which has led to the use of fixed-base 
testing as a common approach used to verify the constrained math models. 
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There are two different attachment configurations that may be used, depending 
on the size of the payload. Larger payloads are supported by four lateral 
trunions and a keel trunion. Figure 2.1 illustrates the payload-orbiter boundary 
conditions for larger payloads. Shorter payloads are supported by two lateral 
trunions and a single keel trunion. Ideally, the trunions are restrained in only 
one or two translational degrees of freedom with all other degrees of freedom 
free. 



Figure 2.1: Payload-Orbiter Boundary Conditions 

For payload math models, the Craig-Bampton form is often used. In 
order to experimentally verify a Craig-Bampton substructure model, the com- 
ponent fixed-interface modes must be validated. Measuring the fixed-interface 
modes involves attaching the substructure to a test-stand in a manner such that 
the appropriate translational displacements at interface degrees of freedom are 
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zero. Therefore, a large and rigid test-stand has to be constructed and tested to 
ensure that its fundamental frequency is above the frequency range of interest 
of the substructure. This requirement is to ensure that there is no coupling 
between the test stand and the test article. Construction of such a test-stand is 
usually very expensive. Frequently, the test stand will not be rigid enough and 
test-stand dynamics must be accounted for by performing additional modal 
tests or finite element analysis of the test-stand. The terms related to the 
constraint modes in the Craig-Bampton model are seldom test-verified. 

Additional difficulties involved with fixed-base testing make the ap- 
proach very difficult and impractical in many cases. Reference [10] describes 
some of the difficulties encountered in the fixed-base test of a Space Station 
module prototype; mainly that the interface degrees of freedom could not be 
properly represented, which introduced nonlineanties into the system at the in- 
terfaces. Alternate test methods have been suggested using free-free modal test 
data to verify the constrained modes. Since a rigid test-stand is not required, 
free-free modal tests are generally less complicated and cheaper to conduct than 
a fixed-base modal survey. However, the test is performed in a different con- 
straint condition than that used for the coupled-loads analysis, so additional 
information must be measured at the constraint degrees of freedom, since free- 
free modes alone are not sufficient to validate a constrained model. 

A mass-additive technique was used as an alternative to fixed-base 
testing to derive constrained modes from free-free modes to verify a Space 
Station prototype module and a Space Shuttle payload [11]. Large masses were 
added at the interface locations and a superposition of the mass-loaded normal 
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modes was used to obtain the constrained modes. The additional masses at the 
interface locations lowered the frequency of the local “trunion modes” into the 
bandwidth of the global modes of the structure. The additional masses also 
allowed the interfaces to be exercised more than during a free-free test. This 
method was found to work better for structures having well-spaced natural 
frequencies than for those with high modal densities, like the Space Station 
module. Drawbacks of this method are that a large number of mass-additive 
modes are required to derive the fixed-base modes, which increases the difficulty 
of the correlation task, and that it may be difficult to determine the proper size 
of the interface masses. Reference [12] contains details of the free-free modal 
test that was used to obtain the experimental data and also a description of 
the problems encountered during the fixed-base modal test of the Space Station 
module. 

The residual flexibility approach has also been used as a technique 
for deriving constrained structural modes from free-free modes supplemented 
with residual dynamics information. A successful application of the residual 
flexibility approach to test-validate a Space Shuttle payload is given in Ref- 
erences [13, 14]. The residual flexibility information was obtained from drive- 
point frequency response functions at each of the payload-shuttle interface loca- 
tions by subtracting the synthesized response from the measured response. The 
payload FEM was then adjusted to match the free-free modes and the residual 
flexibility terms. Admire, et al., also used the residual flexibility approach to 
verify math models of the Space Station prototype module and a Space Shut- 
tle payload [15]. The same approach utilized in Ref. [13] was followed, and 
it was shown that the residual flexibility method worked well and was bet- 
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ter suited than the mass-additive approach for the Space Station module. It 
was concluded that the approach can be used as an alternative to constrained 
boundary testing, but measurement techniques need to be improved to better 
measure the required frequency response functions associated with the residual 
information. 

A technique called Interface Verification Testing has also been used 
as an alternative to fixed-base testing for verification of Space Shuttle payload 
math models [10, 16, 17]. Free-free data is supplemented with additional data 
obtained by exciting a single interface degree of freedom while the structure is 
is freely suspended. The additional information that is obtained in this test 
provides the residual flexibility and is reflected in the location of the zeros of the 
inertance frequency response function. This set of data can be used to verify 
the constrained payload math model. The authors also mention the problem of 
accurately measuring the residual information, since the residual information 
is typically on the order of 1.0 E ® in. /lb, and all of the modal curve fit error 
is contained in the residual curve fit. 

Recent attempts have also been made to design interface connections 
that simulate the payload attachment mechanism. The concept is to have 
the test fixture, or flexure, designed to be very stiff in the restrained trunion 
degrees of freedom and flexible in the un-restrained degrees of freedom. Chung, 
Semaker, and Peebles [18] suggest that to ensure that the restrained and un- 
restrained degrees of freedom do not couple, the axial stiffness should be at least 
three orders of magnitude greater than the stiffness in bending directions. The 
flexure should also be designed to transmit the load properly and to minimize 
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interaction with the test article. Miihlbauer, Troidl, and Dillinger have designed 
a flexure using cylindrical rods with flexible sections at each end that act as 
elastic hinges for an interface degree of freedom [19]. One such flexure is used at 
each trunion interface attachment and these flexures have been used for several 
verification tests with good results reported. 

2.2 Structural System Identification 

Structural system identification is the process of using experimentally 
acquired data to obtain some sort of useful model of the test article that char- 
acterizes the dynamics of the structure. System identification algorithms have 
been formulated in both the time domain and the frequency domain, and each 
has its own advantages and disadvantages [20]. A comparison and evaluation 
of several identification algorithms and the test requirements needed to sat- 
isfy the assumptions of each as applied to the Galileo spacecraft is given in 
Reference [21]. 

The most important advantage of time-domain techniques is that 
they are generally better conditioned numerically than an equivalent frequency- 
domain implementation [22]. This is believed to be due to the power to which 
the frequency values are raised in the frequency-domain equations of motion, 
which increases the dynamic range of the numerical data. For this reason, 
most frequency domain estimation methods generally work better for narrow 
frequency band analysis. Frequency-domain algorithms based on orthogonal 
polynomials have attempted to alleviate this problem [23]. Frequency-domain 
algorithms have the advantage that they can compensate for the effects of out- 
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of-band modes, or residual dynamic effects, while time-domain methods cannot 
represent residual effects [24]. 

Experimental models have traditionally been obtained by experimen- 
tal modal analysis techniques, which attempt to identify the modal parameters 
of the structure, that is, the natural frequencies, damping factors, natural mode 
shapes, modal participation factors, and the residual mass and stiffness [5]. The 
response of the system is described in terms of a linear superposition of the 
characteristic solution of the differential equations. The modal parameters are 
then obtained by curve-fitting the measured data. 

Another approach to obtain the experimental model is to estimate the 
system matrices, that is, to identify the mass, damping, and stiffness matrices 
directly from measured system responses. This approach is referred to as Direct 
System Parameter Identification (DSPI) [25]. The matrices are estimated di- 
rectly from the measured input and output data. An eigenvalue decomposition 
of the estimated matrices results in the desired modal parameters. 

2.2.1 Time-Domain Identification Algorithms 

For time-domain system identification, much effort has recently been 
put forth to develop state-space realization procedures. Most notable among 
these is the Eigensystem Realization Algorithm (ERA) by Juang and Pappa [26] 
and its derivatives [27, 28]. A widely used time-domain identification algorithm 
in the modal test community is the Polyreference algorithm [29]. In Ref. [30], 
Juang relates the algorithmic properties of the Polyreference algorithm to a 
state-space canonical-form realization. 
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Since all vibration test data is discrete in nature, all time-domain 
linear system realization procedures begin with the assumption that a finite 
order discrete state-space model of the system exists, and that it is of the form 

x(k + 1) = Ax(k ) 4- Bu{k) 

y(k ) = Cx(k) + Du(k) (2.2) 

where x is the N x 1 state vector, u is the m x 1 input vector, y is the l x 1 
output vector, A is the N x N state transition matrix, B is the N x m input 
influence matrix, C is the lx N output influence matrix, D is the Ixm matrix 
corresponding to direct input/output feedthrough, and k is the time sample 
index. Note that the C and D matrices in Eq. 2.2 are not the same as the ones 
in Equation 2.1. The solution for the output y{k ) to an arbitrary input u(k) is 
given by 

y(k) = Y.y(k- <)«(■) (2.3) 

i=0 

where Y ( i ) are the discrete-time impulse response functions, which are referred 
to as Markov parameters. The Markov parameters are related to the state-space 
matrices by 

" { CA'-'B t > 0 < 2 ' 4 > 

Given a measurement sequence of Markov parameters, the system realization 
problem is to find a minimal-order realization of the state-space matrices that 
best approximates the given Markov sequence. 

The main problems in system realization are determination of the 
model order and the state-space parameters. The Eigensystem Realization 
Algorithm provides a systematic approach for determining a minimal-order 
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model of & given accuracy and for determination of the discrete state-space 
model. The algorithm uses a discrete-time shift of the Markov parameters to 
form a Hankel matrix, which is defined as 


H qd (k) = 


Y(k + 1) 
Y{k + 2) 


Y (k + 2) 
Y{k + 3) 


Y{k + q ) Y(k + q + l) 


Y(k + d) 

Y(k + d + 1) 

Y(k + q + d - 1) . 


(2.5) 


where q and d are arbitrary integers. The rank of Hgd( 0) is estimated from a 
singular value decomposition, and the largest N singular values are selected to 
determine the model order. The discrete state-space matrices are formed from 
the left and right singular vectors and the Hqd( 1) Hankel matrix. Once the 
state-space matrices are known, the modal parameters can be extracted. 


Su and Juang [31] developed a time-domain algorithm for system 
identification at the substructure level. A procedure is presented to determine 
and assemble substructure Markov parameters. Using the Markov parame- 
ters, substructure transfer functions can be computed and used to determine 
substructure state-space models. A procedure is also outlined to couple the 
substructure state-space models to obtain an analysis model for the entire as- 
sembled structure. The resulting model describes the dynamics of the substruc- 
tures when the appropriate interface compatibility and equilibrium conditions 
are enforced. The authors note that to enforce compatibility and equilibrium 
conditions at the substructure interface locations, co-located sensors and actu- 
ators are needed at all of the interface degrees of freedom. 


The common basis normalized structural identification (CBSI) pro- 
cedure developed by Alvin and Park [32] uses the Eigensystem Realization 
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Algorithm, or any other state-space-based system identification technique, as 
a precursor to determine a second-order structural dynamics model from a 
minimal-order discrete state-space model. The authors point out that for the 
solution of the system realization problem, there are an infinite number of 
equivalent realizations for the given data. A transformation is developed to de- 
termine the underlying physically based structural model from the mathemati- 
cal state-space identified model. The realized model, after the transformation, 
has a one-to-one correspondence with the physical parameters of the system 
and can be used to determine the mass, damping, and stiffness matrices of the 
system. 


2.2.2 Frequency-Domain Identification Algorithms 


Leuridan, et al., present a frequency-domain DPSI approach that uti- 
lizes a straightforward linear estimation to identify the system matrices [25]. 
Beginning with the frequency-domain equations of motion 


' M C K] 


-u 2 X(u) 

juX(u) 

*(«) 


} = DF{ u) 


( 2 . 6 ) 


the system matrices are solved for directly. However in its present form, Eq. 2.6 
cannot be solved very easily. A technique developed in Ref. [33] is used to 
rewrite the matrix of unknowns as a vector and cast Eq. 2.6 in the more familiar 
and readily solvable form Ax = b. However, the resulting equations are highly 
ill-conditioned, and numerical techniques should be used to minimize the effects 
of the ill-conditioning. An advantage of the method is that it uses all of the 
measurement data from multiple inputs simultaneously in the identification 
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procedure. By doing so, it takes advantage of redundant information in the 
data and reduces the influence of measurement noise errors. However, the 
authors point out that if the input frequency data does not fit the model or 
does not contain enough information to form a complete model, the method 
will become unstable and diverge. 

One of the earliest algorithms to identify the system matrioes relates 
the matrices to measured complex modes and frequencies [34]. A similar ap- 
proach is also presented by Balmes [35]. The equations of motion must be 
cast in first-order form to include the effects of damping in the ensuing eigen- 
problem. Just as in the second-order eigenvalue problem, two orthogonality 
conditions result. After some algebraic manipulation and employing the two 
orthogonality conditions, the following expressions are arrived at for the mass, 
damping, and stiffness matrices: 

M = (tf A V T )-\C = -MtyA 2 V T M, and K = -(^A“ 1 ^ T )" 1 (2.7) 

where 'P is the normalized displacement partition of the complex modal ma- 
trix, and A is a diagonal matrix of the complex eigenvalues. An expression, 

= 0, referred to as the propemess condition, is also presented in the 
derivation of Equation 2.7. The propemess condition implies that the mea- 
sured mode shapes are consistent with measured displacement quantities, and 
this constraint must be satisfied in order for the solution of the inverse prob- 
lem to be physically consistent in the determination of the mass and stiffness 
matrices. The algorithm uses the pole and residue information obtained from 
some other algorithm to estimate the complex mode shapes and frequencies. 
The scaling of the complex modes is accomplished through information associ- 
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ated with a drive-point transfer function. A least-squares estimation is used to 
minimize the phase of the complex mode shapes since it is assumed that the 
complex modes should be mostly real. To enforce the propemess condition, a 
constraint minimization problem must be solved for the final estimation of the 
complex modes before the system matrices are computed. In the current for- 
mulation, this algorithm allows for only a single input to be used in the scaling 
procedure. However, most vibration tests employ multiple inputs, so some sort 
of yet-to-be-defined averaging algorithm of the co-located sensor information 
is required before this method can be used with multiple-input test data. 

A frequency-domain direct identification algorithm developed by Lem- 
bregts [36, 37] uses a state-space formulation to solve for system matrices. The 
equations of motion are written in first order form as 



where A x = M~ 1 C, A 2 = M~ l K, and D is the force distribution matrix. A 
Laplace transformation is used to derive an expression to determine A\ and A 2 
from the measured frequency response data. The ensuing first-order algebraic 
eigenvalue problem 



is then solved to identify the modal parameters. The algorithm does allow for 
simultaneous use of multiple input data and also can include terms to compen- 
sate for residual effects. Before the identification procedure begins, an optimal 
model order is selected and a singular value decomposition is performed on 
the measured frequency response data to condense the data to the appropriate 
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model order. In order to compensate for noise in the data, measurement error, 
and residual effects, the model order is selected to be greater than the number 
of physical modes present in the data . In this algorithm, an assumption is 
made that the number of observable modes is less than or equal to the num- 
ber of output measurement locations. This is a fundamental limitation of this 
technique. In test situations, this implies that the frequency band over which 
data is being measured must be narrow enough so that there are fewer modes 
than output measurements. 


Recently, a Unified Matrix Polynomial Approach (UMPA) [38]— [40] 
has been developed to show that many modal parameter estimation methods 
can be reformulated in a consistent framework. The Unified Matrix Polynomial 
Approach is a direct parameter identification method and allows multiple-input 
systems. The model is described in the Laplace domain by a rational fraction 
polynomial of the form 

(d^‘) {*(•)} = (fw) (210) 

where [Aj] and [5j] are the unknown UMPA coefficient matrices, n and m are 
the model orders of the input and output, and s is the Laplace variable (s = ju). 
The coefficient matrices are estimated from the measured input-output data 
using a total least squares estimation having the form 


Aq A\ Ai Bq j 


X(u) ■ 
ju)X (u) 

-FM . 


( 2 . 11 ) 


where /1 0 — K, A, — C, A 2 — M. and B 0 = D. A description of the total 
least squares method is given in Chapter 5. Once the coefficient matrices 
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are estimated, the modal parameters are estimated by solving the algebraic 
eigenvalue problem of Equation 2.9. 

Frequently, in structural dynamics testing the number of output mea- 
surements is less than the number of modes in the frequency bandwidth over 
which the test is being performed. As previously mentioned, some identification 
algorithms require that the number of sensors be equal to or greater than the 
number of observable modes. To get around this problem of spatial truncation, 
Alvin, Peterson, and Park [41] developed a method to determine minimal-order 
mass and stiffness matrices when the number of observable modes is greater 
than the number of available sensors. The resulting minimal-order matrices 
express contributions of all modes that are observable from the available sen- 
sors. The measured physical degrees of freedom are augmented with a set of 
generalized coordinates which possess information about the residual dynamics 
of the system. A singular value decomposition of the residual dynamic matrix 
is used to determine the rank deficiency of the reduced stiffness, and the re- 
quired minimal-rank augmentation is the number of nonzero singular values. 
The minimal-order mass and stiffness matrices are formed by augmenting the 
mode shapes with the singular vectors. A drawback of this method is that it 
assumes that a set of mass-normalized mode shapes can be measured, and it 
is these real normal modes that are used to compute the minimal-order mass 
and stiffness matrices. 

A procedure for identification of modal parameters using both time- 
and frequency-domain procedures is described in Reference [24]. The model 
order is selected using time-domain procedures, and the modal parameters axe 
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estimated by using frequency-domain procedures that are able to represent 
residual effects. To determine the number of modes present in the data, the 
singular values of the frequency response function matrix are examined and a 
modal indicator function is also used. The size of the Hankel matrix is selected 
by observing the convergence of its singular values while varying the number of 
rows and columns of the Hankel matrix. Once the singular values converge and 
the size of the Hankel matrix is chosen, the model order is varied from 1/2 to 4 
times the number of modes observed in the data. Modal quality indicators are 
then used to eliminate poles and to study the convergence of modal frequencies 
and damping ratios. The modal vectors are extracted from the converged poles 
in a least squares frequency-domain curve fit. The frequency-domain curve fit 
includes an improper term as an unknown, which compensates for the fact 
that discrete state-space models cannot represent the effect of modes above the 
Nyquist frequency. The resulting minimal-order model achieves good accuracy 
in both the time domain and the frequency domain. 

2.3 Conclusions 

This chapter has presented numerous testing methods and structural 
system identification algorithms. Each one has its own inherent advantages and 
disadvantages. The selection of one method or algorithm over another depends 
in large part on the structure under investigation, the purpose for which the 
identification is being performed, and on the experience of the structural dy- 
namidst. There is an obvious need for a new identification procedure that can 
overcome many of the previously discussed shortcomings. One such algorithm 
is discussed in the next chapter. 



Chapter 3 


A Proposed Frequency-Domain Substructure System 
Identification Algorithm 


The Substructure System Identification (SSID) algorithm is a new 
frequency-domain identification method that can be used to obtain a linear 
viscous-damped model of a substructure 1 . Every interface degree of freedom is 
either actively excited by a shaker with the input measured or is supported by 
the test stand with the reaction forces measured. The substructure is “isolated” 
from the test stand by measuring the reaction force; thus the substructure can 
be identified. In addition, test stand dynamics are automatically taken into 
account, so there is no need for a separate modal test or finite element model 
of the test stand. This procedure makes it possible to obtain not only the 
fixed-interface modal data for a Craig-Bampton substructure model, but also 
the data associated with the constraint modes. 

An overview of the theoretical derivation of the proposed SSID algo- 
rithm, which is a two-step identification process, is provided in this Chapter. 
A complete description of the algorithm is given in Reference [4]. 


'Although designed especially to facilitate accurate testing of substructures, the algorithm 
can be considered to be a “general-purpose” structural system identification algorithm. 
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3.1 Identification of M 1 C and M l K 

Assume that the substructure has viscous damping and that the total 
number of motion transducers (accelerometers) is at least twice the expected 
number of normal modes in the frequency range of interest 2 . Every interface 
degree of freedom is to have a co-located actuator/sensor (force/accelerometer) 
pair. In addition, there are to be motion sensors (accelerometers) at selected 
interior degrees of freedom. 

Let the equations of motion in physical coordinates and the output 
equation be 


Mx + Cx + Kx = Dp(t) 

y = Sx 


(3.1) 


where x € R Nz is the displacement vector; p e R Np is the input force vector; 
y g R n v is the output measurement vector; M, C, and K are the system mass, 
damping, and stiffness matrices; D is the force distribution matrix, and S is 
the acceleration sensor distribution matrix. For the present discussion, we will 
assume that the above A x -degree-of-freedom model represents a reduced-order 
model of the structure. Let the coordinates be partitioned in the following 
manner: 


/ 



(3.2) 


2 This restriction may be relaxed by the method described in Section 6.4. 
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where / stands for forced DOFs (i.e., DOFs where there is an active force 
input); r stands for reaction DOFs (i.e., interface DOFs where the tested sub- 
structure reacts against the support structure); i stands for interior DOFs 
(i.e., not a DOF where and active force is applied or a reaction is measured); 
and b stands for boundary DOFs, the combination of /-coordinates and r- 
coordinates. These sets of coordinates are illustrated in Fig. 3.1. 



Figure 3.1: Substructure Model - Vibration Test Configuration 


Equation 3.1a can be written in the following partitioned form (damp- 
ing is omitted here): 
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Let us consider the complex frequency response of the substructure 
due to excitation at frequency a;*, but with the interior DOFs force-free. Then, 




(3.4) 
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(From Eq. 3.4 onward, the vectors can be complex.) The complex displacement 
response can be written as 


x k (t) = (3.5) 

and similar expressions can be obtained for velocity, etc. Then, the frequency- 
domain version of Eq. 3.1a can be written as 

M + C + (^h) *] AM = [ D ' Dr 1 { RM } (3,6) 

The experimental data input to the algorithm is complex, but the 
system matrices to be identified are real. To insure that real matrices will be 
determined, a procedure from Ref. [42] is used. After some algebraic manipu- 
lation, the following equation is obtained for estimating the matrices C, K, Dj 
and D r : 


[C K D f D r ] 


ft 


ft 


Hof 

j<jJ 

H af 


L-W 

-Wff] 


9 

9 


Hal 
j U> 
H af 


L— W 2 J 


*[H„] 

-U[Hrf) -mrf) 


= [ ] (3.7) 


where 


C = M~ X C , K = M~ l K , Df = M l D f , and D r = M l D r (3.8) 
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and the H matrices are measured frequency response function (FRF) matrices 
and &[•] and 9f[*] denote the real part and the imaginary part of the given 
quantity, respectively. The data used in Eq. 3.7 is stacked in the following 
manner: 

J37^°i/i( c *''i) "■ (^i) ■“ ju NiM> Ha\f N f (wn u ) 

••• ^H a2 f Nf (u i) ••• j^-H a7 f N{ (u Nu ) 

~jb[Ha Nx f Nj (^l) *•* 

(3.9) 

The other partitions of the data in Eq. 3.7 are formed similarly. 

A least squares solution or a total least squares solution of Eq. 3.7 
is required. A further discussion of least squares solution procedures is given 
in Chapter 5. The effects of noise on the input and output signals can be 
minimized by averaging the FRFs in Eq. 3.7 rather than the force and response 
spectra of Eq. 3.6. 

3.2 Identification of M, C, and K 

In the previous section, an algorithm was described for identifying the 
system matrices C, K, and D, which are defined by Eqs. 3.7. From the system 
matrices C , K, and D, we wish to determine the system matrices Af, C, K, 
and D, especially the first three. The first step is to perform an eigensolution 
using the matrices identified from Eq. 3.7. Let N s = 2 N x , and let us define the 
state variable 
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z = 


( 


X 

X 


) 


and state matrices 



C I 
I 0 


K 
0 -I 



(3.10) 


(3.11) 


Then the following eigenproblem is solved for the complex eigenvalues -V and 
the complex eigenvectors 6 r - 


[\ r A s + B s ]6r = 0 r = 1, ... ,7V, (3.12) 


To determine the system matrices M, C, and K, a mode-superposition 
representation of the complex frequency response can be employed. Using the 
eigenvectors, 6 r , let the mode-superposition solution for the states z be 

z(t) = E^7rW = er(t) (3.13) 

r=l 

where orthogonality holds in the following form: 


0 T v4,0 = diag(or) , 0 t jB, 0 = diag(fcr) (3- 14) 


with 

K 0 
0 -M 

The following modal-response equations are obtained: 


A,= 


C M 
M 0 


, B.= 


(3.15) 


2r7r(t) + = ^£*.P(0 


(3.16) 
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The mode-superposition solution leads, eventually, to the following acceleration 
frequency response: 


A = -<4x t = E (r) K {*•}, WH op ‘ ( 317 > 

r=l ya r / L J 


where 


Pk 

R h 


= r fm \ 

\ R(u k ) j 

\ju k - K) 


(3.18) 

(3.19) 


As in deriving Eq. 3.7, assume that averaged frequency-response data 
are available at N u frequencies, so Eq. 3.17, averaged at each of these frequen- 
cies, leads to the following equation: 

N. 


G af - (-J Rkr {^x} r {4} r \DfD r ] q 


II 

I 

r / 


(3.20) 


If Eq. 3.20 is postmultiplied by Cl} , the following expression for the accelerance 
FRF is obtained: 

N, 


H',M - E (r) K W, {Mn [d,d ' ] H 

r= i V u r / L J 1 


II 

r 

r I 


(3.21) 


Now, the left-hand and right-hand sides of Eq. 3.21 are matrices of dimension 
N x x Nf. Let fj indicate the jth column of each side of this equation. Then, 
write the j>th column of Eq. 3.21 in the form 

- 1 (£) K (M, W!1 { \"; h } < 322 ) 

Finally, let Eq. 3.22 be repeated for each of the N u frequencies and Nf forces, 
and the resulting equations stacked vertically to form the following equation: 
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En E12 ••• Ei t N, 

E Nui ,i E Nu , 2 . . . E Nui>n , J h ( 1/ai 

\ *• 

Eu E12 ... Ei'N. 1 /a.N, 

En„,\ En u< 2 • • ■ En„,N. J J 

(3.23) 

where 

[Ek r]f, = [** {4} r {4}^] {D fj + D r H rfj } (3.24) 

Equation 3.23 is the key equation that is required to estimate the 
system matrices A/, C, and K. It is used to obtain least squares estimates of 
the N a modal parameters a r . The corresponding modal parameters b r can then 
be computed from 

b T = — A r Or r = 1 , • . • , N s (3.25) 

With these values of dr and b r , the state matrices A s and B„, defined by 
Eqs. 3.15, are computed by using Eqs. 3.14, written in the form 

A, = 0 _T diag(ar)0 _1 
B„ = 0 -T diag(6 r )© -1 



(3.26) 
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Finally, the system matrices M , C, and K are obtained by referring to Eqs. 3.15 
and extracting the appropriate partitions of the A a and B t matrices that are 
obtained from Eq. 3.26. 



Chapter 4 


Reduced-Order Analytical Models 


As was discussed in the previous chapter, the output of the Substruc- 
ture System Identification algorithm is the set of characteristic matrices that 
represent the substructure. Since data can be measured only at a limited num- 
ber of locations, it is not possible to measure all of the degrees of freedom of 
a continuous structure. The resulting model is spatially incomplete and can 
be considered to be a “reduced-order” model of the infinite degree of freedom 
structure. To relate the resulting SSID matrioes to some analytical reduced- 
order models, several analytical model-reduction methods are discussed and 
compared in this chapter. 

Model reduction is an attempt to reduce the size of an analytical 
model but still retain the essential dynamic characteristics of the model. An 
overview of model reduction methods is presented in Section 4.1. In Section 
4.2, the Guyan reduction technique is discussed. Section 4.3 describes various 
other reduction methods in use. Finally, a Craig-Bampton type reduced-order 
model that is suitable for test-analysis correlation is presented in Section 4.4. 

4.1 Review of Model Reduction Literature 

The finite element method has become a very useful tool in the anal- 
ysis of modem structures. Stress analysis, failure analysis, internal loads pre- 
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diction, and design modification usually require very fine discretizations, which 
can lead to finite element models on the order of tens of thousands (or more) 
degrees of freedom. Coarse finite element discretizations are not suitable, since 
they usually do not possess sufficient detail for accurate representations of the 
mass and stiffness of the structure. 

To verify that the FEM is sufficiently accurate to predict the struc- 
ture’s response, the FEM must be test validated. The accuracy of the FEM 
is often assessed by comparing the modal parameters of the analytical model 
to those extracted using vibration test data. Test and analysis natural fre- 
quencies can be compared directly. However, the modal vectors cannot be 
readily compared since the FEM will have many more degrees of freedom than 
the test configuration will have accelerometers. In order to compare the FEM 
results with the test results, the test vectors need to be expanded to the an- 
alytical space or the analytical vectors need to be reduced to the test degrees 
of freedom. The former has a severe drawback in that errors present in the 
FEM model that is to be validated are introduced into the test modes that 
are assumed to reflect the true structure. This corruption of the test data will 
generally lead to errors in the model validation process. In addition, there are 
greater computational costs associated with model expansion than with model 
reduction. 

The second approach results in a reduced representation of the FEM, 
a test-analysis model (TAM). This leads to a one-to-one relationship between 
the TAM degrees of freedom and the number of accelerometers in the test 
configuration. A number of procedures for generating reduced-order models 



32 


have been developed in the past. The reduction methods are based on Ritz- 
type transformations of the form 

B = T T AT (4.1) 

where A is the original matrix, B is the new matrix, and T is the transfor- 
mation matrix. The differences between methods lie in the transformation, or 
interpolation shapes, used to represent the motion of the non-instrumented, 
or omitted, degrees of freedom in terms of the measured, or retained, degrees 
of freedom. The measured degrees of freedom will be referred to as the a-set 
(active degrees of freedom) and the omitted degrees of freedom will be referred 
to as the o-set. 

The model reduction can be done at the system level or at the com- 
ponent level using component mode synthesis (CMS) methods. The two ap- 
proaches used most often in system model reduction are Guyan, or static, re- 
duction and modal reduction. The modal reduction approach lends itself very 
easily to component-mode-synthesis based substructunng methods. Static and 
modal reduction methods differ in both accuracy and robustness, terms which 
axe defined in Reference [43]. Reference [43] defines accuracy as the measure of 
the capacity of the reduced-order model to predict the modal frequencies and 
mode shapes of the FEM. Robustness is defined as a measure of the TAM’s 
ability to show orthogonality of test modes when the FEM has inaccuracies. 

4.2 Guyan Reduction 

The simplest and most straightforward reduction procedure is the 
Guyan, or static, reduction method [44]. This method is very useful for gener- 
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ating test-analysis models since the measured degrees of freedom can be selected 
as the ones to be retained in the reduction process. However, the Guyan re- 
duction method is quite sensitive to the selection of omitted degrees of freedom 
and often results in poor accuracy if there is inertia associated with the omitted 
degrees of freedom. Various authors have devised ways for automated selection 
of the best degrees of freedom to retain [45, 46]. 

The principal assumption of the Guyan reduction procedure is that 
inertial forces associated with omitted degrees of freedom are negligible in com- 
parison with the elastic forces transmitted to the omitted degrees of freedom 
by the motion of the retained degrees of freedom. The mathematical statement 
of the above assumption is based on solving a static problem of the form: 



If there are no loads applied at the omitted degrees of freedom, the upper 
partition of Eq. 4.2 can be solved, giving 

Xo = -K^KoaXa (4.3) 

This yields the following Guyan transformation matrix from the original degrees 
of freedom to the a-set degrees of freedom: 

T static = ~ K °° Koa (4.4) 

The rth column of the transformation matrix represents the static displacement 
of the structure when the rth a-set degree of freedom has a unit displacement 
and the remaining a-set degrees of freedom have zero displacement. This dis- 
placement type of shape is commonly referred to sis a constraint mode. If 
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the modal shapes in the frequency range of interest can be expressed as lin- 
ear combinations of the constraint modes, the Guyan model will be a good 
approximation. 

In the above transformation, it was assumed that there are no forces 
applied to the omitted degrees of freedom. This is not an accurate assumption if 
the omitted degrees of freedom possess significant mass. If the omitted degrees 
of freedom possess large mass-to-stiffness ratios, the omitted inertia terms will 
be nontrivial, and will cause errors in the natural frequencies and mode shapes 
of the reduced-order model. 


4.3 Other Model-Reduction Methods 

Since the introduction of the Guyan model reduction method, various 
attempts have been made to improve upon the Guyan reduction method. A 
method developed by O’Callahan [47], called the Improved Reduced System 
(IRS), includes a static approximation to the dynamic inertia terms of the 
omitted degrees of freedom: 


T[RS — TftaUc "I" T dynamic 


(4.5) 


where 


T dynamic — 


^oo [Moa + Af ooT static] statu: 

I 


(4.6) 


where M^uc = Static = Tj^KT^ and M and K are the 

original FEM mass and stiffness matrices, respectively. The resulting eigen- 
solution of the reduced-order IRS model is generally more accurate than that 
obtained using the original Guyan method, but it is still an approximation 
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of the original FEM dynamics. The original IRS method is applicable for 
global, or system, models and is not directly applicable for substructure mod- 
els. Flanigan describes an extension of the IRS method which works with 
dynamic substructures to create substructure level TAMs [48]. 

An analysis of the robustness of the IRS method has been given by 
Gordis [49]. He points out that even though the IRS TAM may provide better 
estimates of the mode shapes and natural frequencies than the Guyan method, 
the static TAM may provide better orthogonality of the test modes. This is due 
to the approximation of the inertia terms by the IRS method and is dependent 
upon the degrees of freedom that are retained. The retained degrees of freedom 
should be selected so that the lowest frequency of the modes of the omitted 
system is above the frequency range of interest; otherwise the approximation 
of the inertia term will be poorly conditioned. 

Kammer introduced a TAM methodology that is a modal reduction 
method [50]. The advantage of this method is that it produces an increase in 
accuracy of the natural frequencies and mode shapes. This method provides 
an exact transformation between the FEM degrees of freedom and the TAM 
degrees of freedom by using the target modes of interest as the basis of the 
transformation to relate the omitted degrees of freedom to the retained degrees 
of freedom. The modal transformation is given by 

Tmodal = Qt&L (4.7) 

= (*£*<.)■'*£ (4.8) 

where the columns of 4> t are the target modes, and is the generalized inverse 
of the target modes partitioned to the retained degrees of freedom. This same 
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transformation was also presented by O ’Callahan as the System Equivalent Re- 
duction/Expansion Process (SEREP) [51]. The resulting reduced-order model 
is in physical coordinates. When the number of target modes used to form the 
modal transformation is equal to the number of active degrees of freedom, the 
modal mass and stiffness matrices contain an exact description of the FEM’s 
dynamics for the modes of interest, that is, the natural frequencies and mode 
shapes are exactly the same as those of the original model. To ensure that 
Eq. 4.8 is numerically well conditioned, the active degrees of freedom chosen 
should render the target modes linearly independent and observable. 

The static TAM is more robust than the modal TAM in situations 
where the modes of the test article differ significantly from those measured 
during the vibration test. This is true even when the modal TAM accurately 
predicts the natural frequencies. The modal reduction method uses a small 
number of FEM mode shapes to develop the transformation matrix and this is 
more limiting than the static methods which use constraint modes for each a-set 
degree of freedom [43]. Rammer [52], and Bhatia and Allemang [53], point out 
that in some cases, the use of a modal TAM in test-analysis orthogonality and 
cross-orthogonality computation can result in larger off-diagonal terms than 
those produced by a “less accurate” static TAM. The authors hypothesize that 
the discrepancies between the test mode shapes and the FEM mode shapes are 
due to the modal TAM’s poor representation of residual modes. Here, residual 
modes are those modes that are not included in the transformation matrix. 

In an effort to extend the capabilities of the modal TAM method to 
include an improved representation of the residual dynamics, Rammer devel- 
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oped the a new method called the hybrid TAM [52]. The hybrid TAM combines 
the robustness of the static TAM with the accuracy of the modal TAM. Using 
the target modes, an oblique projector matrix is formed to divide the vector- 
space containing the TAM dynamics into two complimentary subspaoes. The 
projector matrix, P, can be expressed as 

P = ** T T^ l MT modal (4.9) 

where 4> is the matrix of FEM mode shapes, is the modal transformation 
matrix (Eq. 4.7), and M is the original FEM mass matrix. Note that the last 
three terms are just the modal TAM mass matrix. The transformation matrix 
for the hybrid TAM can be expressed as a combination of the Guyan and modal 
transformation matrices 


T hybrid — Trtatjc -f- ( Tjnodal Tstati^P (4.10) 

This TAM preserves the exact representation of the target modes and also 
produces a more accurate representation of the residual modes. 

In an approach very similar to Ref. [52], Bhatia and Allemang [53] 
also develop a transformation that preserves the exact FEM target modes and 
includes residual mode information. However, their projector matrix is written 
as 

P=[*ta 0 ] [ 4>to 4Va] _1 (4.11) 

where 4^ is the a-set partition of the target modes and 4> ra is the a-set partition 
of the residual subspace. The residual mode information is obtained from the 
system flexibility matrix and the target modes; thus the residual modes are 
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never computed directly. The residual information is obtained from 

3> r = G - $ t [diap(A?)]^ (4.12) 

where G = K~ l is the system flexibility matrix, 4> t are the target modes, and 
[diag(Xi)} is a diagonal matrix containing the target natural frequencies. 

4.4 Craig-Bampton Reduced-Order Models 

Component mode synthesis (CMS), or substructuring, is an alterna- 
tive approach for model reduction. In the CMS approach, the FEM is first 
divided into several smaller substructures. Each substructure is itself reduced 
to a smaller number of degrees of freedom. The reduced degrees of freedom 
include physical coordinates at substructure interface locations and generalized 
coordinates that represent the component modes of the substructure. Once the 
substructures axe analyzed, the components are coupled together to form the 
full system model. 

The Craig-Bampton method, a variant of the Hurty method [54], is 
the most popular CMS method, and its derivation is described in Reference [55]. 
The transformation in the Craig-Bampton method involves component-level 
fixed-interface normal modes augmented with a set of constraint modes. This 
results in a model that is statically exact at the boundary locations. This 
method is usually more efficient and more accurate than a static reduction 
when the same number of degrees of freedom are retained. 
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4.4.1 Craig- Bampton Model Reduction 


The N physical degrees of freedom of the substructure are first par- 
titioned into two sets, the boundary, or interface, degrees of freedom and the 
interior degrees of freedom: 


x = 



(4.13) 


where Xb are the boundary degrees of freedom and Xi are the interior degrees 
of freedom 1 . Omitting damping, the equations of motion are written as 


' Mu m* ] r Xi \ r Ku K ib ] r ** \ _ r u \ 

. M u Mtb J \ x b j + K h ^ J \ x b J \ f b f 

The Craig-Bampton method uses the Ritz transformation 


(4.14) 


x = TcBiQ 


where 


Tcb 1 = 


0 I 


* c =-K~ i 1 K ib 


and q is a vector of generalized coordinates having the form 


Q = 



(4.15) 


(4.16) 

(4.17) 


(4.18) 


The first column-partition of the Craig-Bampton transformation ma- 
trix contains the fixed-interface normal modes. This describes the motion of 


^though originally developed as a component mode synthesis method, the method can 
be considered to be a general model-reduction method by simply letting i — ► o and b — ► 
a. That is, let “interior 51 coordinates be more general “omitted” DOFs and, likewise, let 
“boundary” coordinates be more general “active” coordinates. 
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the interior degrees of freedom relative to the boundary degrees of freedom in 
terms of the normal modes of the substructure with the boundary degrees of 
freedom fixed. The fixed-interface normal modes are obtained by solving the 
eigenproblem for the interior degrees of freedom 

(Ku — L%Mn)<t> r = 0 (4.19) 

or 

KA = MA A n (4.20) 

where the subscript r denotes the rth fixed-interface normal mode. The second 
column-partition of the transformation matrix contains the constraint modes. 
These are similar to the constraint modes described previously. The rth con- 
straint mode is defined by producing a unit displacement at the rth boundary 
degree of freedom with all other boundary degrees of freedom constrained and 
with all interior degrees of freedom unconstrained. 

The final reduced-order matrices obtained are a combination of phys- 
ical (boundary) coordinates and generalized modal coordinates. These gen- 
eralized coordinates cannot be used directly for test-analysis correlation and 
therefore must be transformed into suitable physical coordinates. 

4.4.2 Craig-Bampton Models in Physical Coordinates 

The transformation from generalized coordinates to physical coordi- 
nates for Craig-Bampton models is described by Huang and Craig [56]. The 
interior physical degrees of freedom are partitioned via 

<*>-{£) 


(4.21) 
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where x r axe the degrees of freedom where measurements will be made during 
the vibration test and x 0 are the omitted degrees of freedom. The equations of 


motion can be rewritten as 
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(4.22) 


Then, Eq. 4.15 can be written as 

= T CBl q = 

0 I 

This leads to the transformation matrix Tcb 3 relating the generalized ooordi- 


(4.23) 


nates to physical coordinates to be used during the test: 


Tcb 2 = 



(4.24) 


A similar approach was developed by Kammer [57], but Huang and Craig point 
out that the contribution of the constraint modes to the transformation matrix 
Tcb? was not included in Rammer’s formulation of the transformation. 


The order of the Craig-Bampton model is determined by how many 
fixed-interface normal modes are used to augment the constraint modes when 
forming the transformation matrices. If the desired size of the reduced-order 
model is k, then the number of fixed-interface normal modes retained is n m = 
k — rib. Typically, only the lowest rim fixed-interface normal modes are kept. 
This results in a Craig-Bampton model that accurately predicts the lowest k 
modes of the substructure. However, if the target modes are not the lowest 
consecutive modes, then the analytical modes from the Craig-Bampton model 
will not necessarily correlate well with the target test modes. 
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In order for the Craig-Bampton model to represent the target modes, 
a proper set of fixed-interface normal modes must be selected. This set con- 
sists of the fixed-interface normal modes that contribute most to the target 
modes. The selection process is accomplished in the same manner that was 
used in Ref. [58] to determine which coordinates to select as the rigid-body 
coordinates to solve a semidefinite eigenvalue problem. Gaussian elimination 
or other factorization methods can be performed on 4> n r with pivot selection 
being used to determine which columns of $ nr contribute most to the target 
modes. Pull column and row pivoting is necessary in the event that one of the 
retained interior degrees of freedom lies along a node line of a fixed-interface 
normal mode. It was observed, however, that this selection technique was de- 
pendent upon the scaling of the fixed-interface normal modes and further study 
is needed to remove the scaling dependence. 

A Craig-Bampton test-analysis model suitable for correlation with 
the test modes can be obtained by choosing the fixed-interface normal modes 
that contribute most to the test modes. The final form of the transformations 
defining the Craig-Bampton TAM are given by 

Mcb = T t cb ?Z b JAT C bJcb> (4-25) 

Kcb = '^'cBi^'cBi KTcbi Tcb, (4-26) 

Icb — ^cb^cbJ (4-27) 

where the transformation matrices contain the proper set of fixed-interface 


normal modes. 



Chapter 5 


Computational Aspects of the SSID Algorithm 


This chapter describes some of the computational considerations as- 
sociated with the SSID algorithm. Section 5.1 contains a description of various 
methods of solution for a set of linear, algebraic equations. The ill-oonditioning 
associated with frequency-domain identification procedures and a possible rem- 
edy for this problem are discussed in Section 5.2. The topic of model order 
determination is addressed in Section 5.3. 

5.1 Least Squares Methods 

The solution of systems of linear, algebraic equations is involved in the 
key mathematical steps in the Substructure System Identification algorithm for 
accurately identifying the substructure. The method of solution can drastically 
affect the resulting identified system. The methods of solution used in the 
numerical simulations in this report are the least squares (LS) method and the 
total least squares (TLS) method. A brief explanation of the methods and of 
their differences is given here. 

The least squares method and the more recent total least squares 
method are both mathematical modeling procedures used to solve an under- 
or over-determined linear system of equations of the form 

AX = B (5.1) 
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where A is an m x n data matrix, X is an n x d matrix of unknowns, and 
B is an m, x d matrix of observations. A set of under-determined equations 
results when there are more unknown parameters than equations; for the over- 
determined system, there are more equations than unknowns. Unless B belongs 
to the range of A, R(A) (i.e., the set is compatible), the over-determined system 
of equations has no unique solution. This case, which results when the data 
contains noise, will be denoted by 

AX « B (5.2) 

The approximation symbol is used to emphasize that the data (i.e., A and/or 
B ) may be contaminated by noise. If there is no noise in the data, then the 
equality is valid. The least squares method and the total least squares method 
seek a solution that minimizes the error, e, between the true system model and 
the measured data. The solution obtained depends on the error model and on 
the weighting of the data used by each method. If the noise does not match the 
error model, a biased estimate of the solution will result. An unbiased estimate 
is one that, on the average, neither tends to over-estimate nor under-estimate 
the true solution. 

Figure 5.1 geometrically illustrates the error models used in the two 
methods. By depicting the least squares and total least squares methods as 
measures of goodness-of-fit, the least squares problem is one of minimizing 
the siim of the squared distances along a single coordinate axis. Whereas the 
total least squares problem seeks to minimize the sum of the squares of the 
perpendicular distances of the observed data from the fitted line to provide the 


best fit. 
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Figure 5.1: LS and TLS Error Models 


5.1.1 The Least Squares Method 

Perhaps the best known method of solution for an over-determined 
system of linear equations is the least squares method. For simplicity, consider 
the case when d, = 1. The accuracy of the least squares solution is independent 
of the number of right-hand-side vectors used in the computation. In the 
classical least squares approach, the measurements of the variables in the data 
matrix A are assumed to be free of error; hence all errors are confined to the 
observation vector b. A least squares estimation could also be made assuming 
that the observation vector is exactly known and the data matrix contains 
errors. The least squares method will result in an unbiased estimate if the 
error model of the given data is of either of the following two forms: 

Ax = {bo + A6} (5.3) 

[Aq + A Ai\x = b (5.4) 

where Aq is the true data matrix, bo is the true observation vector, and A A and 
A b are the corresponding errors. In most situations, Eq. 5.3 is the form that 
is assumed, so it will be used here in the following discussion. It is assumed 
that the errors A bi and A Ay, i = 1, . . . , m and j = 1, ,n are un correlated 
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random variables with zero mean and equal variance. This assumes that any 
necessary preprocessing measures on the data such as scaling and centering 
have been performed in advance. 

The classical least squares estimate is equivalent to minimizing the 
sum of the squares of the differences, or the 2-norm, between the measured 
observation vector b and an estimated observation vector b. The least squares 
problem seeks to min \\b - 6|| 2 where b is the orthogonal projection of b onto 
the i?(A). This amounts to perturbing the observation vector 6 by a minimum 
amount Ab so that b = b - Ab can be predicted by the columns of A, that 
is, f> lies in the R(A). The minimum perturbation, Ah, is the called the least 
squares correction. The resulting least squares solution, assuming that A is of 
rank n, is 

x = (A T A)- 1 A T b (5.5) 

In addition to the above normal-equation solution, the least squares solution 
can also formulated using other techniques, such as the singular value decom- 
position or the Q-R decomposition. 

The main assumption made in the classical least squares estimation 
is that errors only occur in the observation vector b and that the data matrix 
A is independent of noise. However, this assumption is frequently unrealistic: 
sampling errors, modeling errors, instrument errors and human errors may 
imply inaccuracies of the data matrix A as well. This will cause the least 
squares solution to yield a biased estimate. 
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5.1.2 The Total Least Squares Method 

The total least squares method [61] was developed to provide esti- 
mates that are generated from a system of linear equations where it is assumed 
that both the data matrix A and the observation matrix B contain noise. The 
error model for the total least squares method is of the form 

[Ao + AA){X} « [B 0 + A B] (5.6) 

The error criterion again assumes that the errors are uncorrelated random 
variables with zero mean and equal variance. 

The total least squares problem is formulated by rewriting Eq. 5.1 as 
a homogeneous system of linear equations 

\A;B\ «0 (5.7) 

The approximation symbol is again used to emphasize that the data is con- 
taminated by noise. If there is no noise in the data, the right hand side 
of Eq. 5.7 will be identically zero. The total least squares problem seeks to 
min || [A] B] — [.A; fi]|| F subject to B € R(A), where A and B are the total 
least squares approximations of A and B required to make the set of equations 
compatible. The Probenius norm, || • || F , is a measure of the size of an m x n 
matrix; similar to the vector 2-norm, which measures the length of a vector. 
Again, for simplicity, assume that d = 1. 

The solution to this homogeneous system of linear equations is or- 
thogonal to the row space of the augmented data matrix, that is, it is equal to 
the nullspace of the augmented data matrix. The nullspace of a matrix consists 
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of all vectors x such that [A; b]x = 0. There are numerous methods that can 
be used to divide the augmented data matrix into a system subspace and a 
null subspace, but the numerically stable singular value decomposition is most 
often used. 

Singular value decomposition methods are based on the following the- 
orem of linear algebra [59]: Any mxn matrix Z can be written as the product 
of an m x m column-orthogonal (or unitary if Z is complex) matrix U, an mxn 
diagonal matrix E with non-negative elements, and the transpose (Hermitian) 
of an n x n orthogonal matrix V , that is, 

Z = UEV t (5-8) 

Assuming that Z is of rank r where (r < min{m,n}), then [60]: 

• The columns of U are called the left singular vectors of Z and are the 
eigenvectors of ZZ T . 

• The first r columns of U form the column space of Z . 

• The last m — r columns of U form the nullspace of Z T . 

• The diagonal elements of E are the singular values (<x<, i = 1, . . . ,n) of 
Z sorted in descending order. The rank of Z is equal to the number of 
non-zero singular values. 

• The singular values squared are the eigenvalues to both ZZ T and Z T Z. 

• The columns of V are called the right singular vectors of Z and are the 
eigenvectors of Z T Z. 
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• The first r columns of V form the row space of Z. 

• The last n — r columns of V form the nullspace of Z. 

Notice that the columns of U and V give orthonormal bases for all four funda- 
mental matrix subspaces. 

The singular value decomposition expresses a matrix as a linear com- 
bination of rank-one matrices. The best rank p(p <r) approximation of matrix 
Z is obtained by retaining the first p terms of this decomposition, which cor- 
responds to the p largest singular values. The error of this approximation is 
the smallest among all rank p matrices, and is equal to the next singular value 
<7p+i when the 2-norm is used to measure the error. 

Let the singular value decomposition of the augmented data matrix 
[A] 6] be given by 

[A]b] = UZV T (5.9) 

where U is an m x m matrix, E is an m x (n + 1) matrix, and V is an (n + 
1) x (n + 1) matrix. If there were no errors in A and 6, b would lie in the space 
defined by the columns of A and the set would be compatible. The rank of 
matrix [A] b] would be n and the corresponding total least squares correction 
would be zero. The total least squares solution could be determined from the 
last right singular vector. However, since A and b contain errors, the matrix 
[>1; 6] is of rank n -I- 1 and there is no nonzero vector in the nullspace of [A] b ], 
N([v4; 6]). To obtain a compatible set, the rank of matrix [A; 6] must be reduced 
to n. The best rank n approximation of [ A ; b] is given by 


[A] 6 ] = utv T 


(5.10) 
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where £=diag(<7i, . . . ,<r„,0). The total least squares minimal correction is 
thus a n +\ = min || [A\b] - [A\b] ||f- The approximate set of equations is now 
compatible and its solution is given by the only vector v n+i that belongs to 
W([A ; S]). The total least squares solution is obtained by scaling v B +i so that 
its last component is -1 or 



-1 

v n +\ 

^n+l,n+l 


(5.11) 


If the rank of the nullspace is greater than one, multiple estimates 
of the solution vector will exist. This occurs whenever a n = <r n + 1 , or more 
generally if a v > Op+x = • • • = °n+i with p < n. This situation also occurs 
whenever the system of equations is under-determined. Instead of accepting 
all the solutions as estimates and then performing some type of averaging, a 
linear combination of the estimates will be used to find the minimum norm 
total least squares solution. 

In the case when v n +i, n +i = 0, the total least squares problem may fail 

to have a solution. This type of problem is referred to as a nongeneric problem 

in the literature. Nongeneric problems occur whenever A is rank-deficient, or 

when the set of equations is highly conflicting. A total least squares solution 

X 

may still be found by imposing an additional constraint ^ -L Un+i* 

A generalization of these results can be made to allow for multiple 
right-hand sides, that is, for d > 1, where B € iT xd . For example, using 
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Eq. 5.7, Eq. 3.7 can be written as 

r (£M(c*) ] 

r , 

[ C Df D r l]l -F(u k ) f « 0 (5.12) 

-R(u k ) 

A{uk) 

The singular value decomposition of matrix [A; B ] can be computed as: 

M B ) = USvT = l U ‘ : [ o‘ l][v 2 \ v“f <513) 

where U is an m x m matrix, E is an m x (n + d) diagonal matrix, and V 
is an (n + d) x (n + d) matrix. When the rank of matrix [ A ; B] is greater 
than n, the set is incompatible. The last d singular values, cr n+1 , . . . , cr n +d, are 
not equal to zeros. The total least squares method forces the equations to be 
compatible by setting the last d singular values to zero. These are the smallest 
singular values and are due mainly to the noise in the data. Again, this is the 
minimal correction to the augmented matrix [A; B]. The right singular vectors 
corresponding to the last d singular values define the nullspace of the augmented 
matrix and form the total least squares solution. After being properly scaled, 
the solution is 

x = -V n Vn' (5.14) 

In practical applications, a threshold value must be set to determine 
which of the singular values should be set to zero. The threshold value depends 
in part on the quality of the measurement data. Estimations for the noise 
level on the data, if such data is available, can be used to determine this 
minimum value. Note that the ratio ^ cannot exceed the precision limit of 
the machine [62]. 
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5.1.3 The Scaled Total Least Squares Method 

In both of the methods previously discussed, it is assumed that any 
necessary preprocessing of the data, such as scaling and centering, have been 
performed in advance. This ensures that the measurement data used in the 
computation is equalized in magnitude and will result in the best unbiased 
solution. 

The data matrix in Eq. 5.12 can be decomposed into two matrices by 
using the total least squares error model given in Eq. 5.6, 

[ X /] H Z« 0 (5.15) 

where 

X H = [ C K Df D r ] (5.16) 

and 

(^)A>(^fc) ) f (j^) E a(^k) 

(i)AiW 

Z = Zo + E =l -FoM > + < -E f (u k ) | (5-17) 

- Ro(vic ) -E r {uj k ) 

. i4(o; fc ) J l E a (u k ) . 

where Z 0 is the true error-free measurement data and E is the measurement 
error matrix. Inspection of Eq. 5.17 reveals that some sort of scaling of the 
measured data is needed before use in the system identification process. For 
the frequency range from 10 Hz to 200 Hz, the multiplicative factors associated 
with the errors range from 1 to 6 x 10 -7 . Therefore, in order to achieve the 
best unbiased estimation, the data matrix Z should be scaled so that the noise 
is row-wise equalized. 
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Tb perform this equalization, Li, et al. [39, 40] utilized a scaling pro- 
cedure that is referred to as the scaled total least squares method (STLS). 
For the examples given in Refs. [39, 40], the STLS method results in better 
estimates than do the least squares method and the total least squares method 
without scaling. However, the scaled total least squares method assumes that 
an estimate of the error matrix is available. In test situations, this information 
is usually not known. A more practical scaling method is thus needed, but 
Eq. 5.17 emphasizes the need for proper scaling of the measurement data. 

It is important to remember that least squares and total least squares 
are only two possible techniques for estimating the unknown parameters of a 
linear system. The total least squares method does, however, give the best 
estimates (in a statistical sense) when the error of the system satisfies the total 
least squares error model [61]. 

5.2 Narrow-Band Data Processi n g 

It has been mentioned previously that frequency-domain identification 
algorithms generally result in better estimations when used as narrow-band 
identification procedures [36], This is believed to be due to the exponent in 
the frequency-domain equations of motion. However, sometimes not all of the 
target modes lie within a narrow frequency band. For instance, the target 
modes may span a frequency range of 2 kHz or more. Terms in the resulting 
equations of motion would vary by a minimum of 8 orders of magnitude, and 
this could lead to the data being ill-conditioning. In the SSID algorithm, the 
first estimation procedure, (i.e., the solution of Eq. 3.7) is more prone to ill- 
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conditioning than is the second step (i.e., the solution of Eq. 3.23) since this 
step uses accelerance, mobility, and receptance frequency response functions. 

In the original formulation of the SSID algorithm [3], all of the mea- 
sured data was to be processed simultaneously. As just pointed out, this could 
lead to numerical difficulties, lb alleviate this problem, a concept referred to 
band processing has been employed. The basic idea behind band processing 
is to work with narrow bands of frequency data individually, and thus avoid 
some of the conditioning problems. In this manner, experimental data span- 
ning several decades could be processed incrementally, for example 500 Hz at 
a time, by the SSID algorithm until all of the data has been included in the 
identification. 

The band processing is implemented in the SSID algorithm as follows. 
First, the experimental data is acquired in the usual manner covering as many 
decades as necessary to obtain all of the desired modes. The data is then 
divided into suitable bandwidths of overlapping frequency bands. The overlap 
allows redundant processing of all of the data; in effect, all of the data is 
processed twice. Then, the SSID algorithm is used to process each frequency 
band individually though the solution of the state-space eigenvalue problem, 
Equation 3.12. From each complex-conjugate pair of eigenvalues, a natural 
frequency and damping factor can be computed by 


II 

1 

P> S3 
— £ 

r = l,...,N x 

(5.18) 

_ 9(A r ) 

Mr — / — 

v'l-e 

rH 

II 

(5.19) 


This step is repeated until all of the frequency bands have been processed. 
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The next step is to separate the structural roots from the compu- 
tational, or noise, roots. In this study, four criteria were used to disting uish 
between the physical modes and computational modes. 

1. To graphically aid in the selection process, the estimated natural fre- 
quencies and the frequency bands are plotted versus the entire frequency 
bandwidth. Any estimated frequencies lying outside of the given fre- 
quency band are discarded. 

2. Roots with negative damping values are discarded. 

3. The modal phase collinearity (MPC) is computed for each complex-conjugate 
pair of mode shapes. Those modes with less than a specified MPC value 
are discarded. The value is dependent upon the level of noise in the 
measurement signal and the quality of the data. 

4. If the root is not repeated in another frequency band, it is discarded. 
Thus, the need for the overlapping frequency bands. True structural 
roots will be estimated in multiple frequency bands, while computational 
ones are highly unlikely to be repeated in multiple frequency bands. 

The flow chart in Figure 5.2 summarizes the steps in the band processing 
procedure. 

The modal phase collinearity is an index developed by Juang and 
Pappa [26] that measures the linear functional relationship between the real and 
imaginary parts of the mode shape coefficients. It based on the assumption that 
for lightly damped structures, the estimated mode shapes from Eq. 3.12 should 
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be close to normal, that is, the phase angle between the complex coefficients 
of the same mode shape should either be 0° or ±180°. The MPC index is 
computed from 



MPC _ faP + gL 

fell 2 + fell 2 

(5.20) 

where 

N x 

(5.21) 


c _ fell 2 - fell 2 
2»3fe 

(5.22) 

and 

a = arctan(|e| + y/l + e 2 ) 

(5.23) 

and || - 

|| represents the Euclidian norm of a vector. An MPC index of near 


unity indicates a normal mode, while a low index indicates a “noise mode.” 
It should be noted that the MPC assumes real normal mode behavior of the 
given structure, which is valid for most structures. 

Once the computational roots are eliminated and the physical roots 
stored, the SSID algorithm can proceed to identify the system matrices just 
as before using the remaining structural roots. Any of the repeated roots and 
corresponding eigenvectors can be selected from a given frequency band for use 
in the second part of the SSID algorithm since the eigenvectors are arbitrarily 
scaled. The scaling used to determine the system matrices arises from the 
solution for the a coefficients in Equation 3.23. 

In the band processing simulations included in the Chapter 6, three 
different approaches where used to estimate the dr values. The first was to 
solve Eq. 3.23 and estimate the d coefficients using all the FRF data in a given 
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frequency band from which the root was selected. This approach is referred 
to as BP1 in Section 6.3. The second approach is similar to the first, except 
that only data near resonance is used in the estimation of the dr values. This 
approach is referred to as BP2. The third approach is to use the data obtained 
near the resonances of all of the selected roots simultaneously in Gq. 3.23; this 
is referred to as BP3. The first two approaches provide a local estimation, 
while the last is a global one since data from multiple frequency bands is used 
simultaneously. 

Band processing also holds the potential of significant computational 
savings. Since the computations are done band by band with a much smaller 
amount of data, each frequency band could be processed simultaneously on a 
parallel-processing computer. For tests where a large number of frequency fines 
are obtained, this could prove to be very cost effective. 

5.3 Model Order Determination 

Perhaps the most important, and most difficult, part of system iden- 
tification is the determination of the proper model order. For structures with 
a very high modal density, there are typically more modes in the frequency 
bandwidth than there are output sensors. For many system identification al- 
gorithms, including the SSID algorithm, this poses a serious problem that they 
are unequipped to handle. In its present formulation, the model order is lim- 
ited to the number of output sensors. This is the only serious drawback of the 
method. 

For the case when there are more modes present than sensors, the size 
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of the model must be expanded. An attempt to expand the size of the identi- 
fication has been developed using pseudo degrees of freedom. The additional 
degrees of freedom are obtained by stacking the data used in Eq. 3.7 differently. 
The input frequency spectrum is divided in half, into a low-frequency spectrum 
and a high-frequency spectrum. In Eq. 3.9, the low-frequency response is then 
stacked on top of the high-frequency response as shown: 


'Ha / 


( W l) * • ■ (“N.) 

jui ^ a N x h (^l) H aNz f N f (cujv. ) 

A ' ‘ ’ (vnJ 


(5.24) 


[ jVfV e+1 ^ a N x f\ l) ^ a N x In J ) J 

where N e = Stacking the data in this manner essentially doubles the order 
of the model. 


Thus, although the lower partition of data is measured at the same 
physical coordinates as the upper partition of lower-frequency data, the system 
modes that contribute to the two sets of data are not identical. The lower 
partition is referred to as data from pseudo degrees of freedom. After Eq. 3.7 is 
solved using the expanded model, the state-space eigenvalue problem, Eq. 3.12, 
is solved. However, the details of the second step of the algorithm, the iden- 
tification of the system matrices, have yet to be resolved. Further research is 
needed in this area. Section 6.4 presents results obtained by processing through 
the first step of the SSID algorithm employing pseudo degrees of freedom. 



Chapter 6 


Numerical Simulations 


This chapter describes the numerical simulations used to study the 
Substructure System Identification algorithm. The simulations presented here 
focus on model-reduction and the SSID algorithm s ability to identify a valid 
reduced-order structural model. This is a problem that should be addressed by 
all system identification algorithms, since spatial truncation is inevitable in test 
situations. A proof-of-concept example and examples using lumped-parameter 
models are given in References [3, 63]. 

In Section 6.1, a description of the analytical model and an overview 
of the simulations are given. The results from several simulations are pre- 
sented in Section 6.2. Results from simulations employing band processing are 
included in Section 6.3. Section 6.4 presents the results when the model or- 
der is expanded by the introduction of pseudo degrees of freedom. The SSID 
reduced-order models are compared to other TAMs in Section 6.5. 


6.1 Model Description and Overview of Simulations 

The analytical model used in the simulations is a 52-DOF FEM of the 
“Payload Simulator” in the Structural Dynamics Laboratory at The University 
of Texas at Austin. This is a model of the structure that was used in a previous 
study to measure the reaction forces at the substructure-test stand interface 
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locations [64]. The physical structure consists of two 60-in .-long aluminum box 
beams connected by two 20-in. cross-beams at either end. 

The FEM consists of 18 nodes and 20 beam elements and is illus- 
trated in Figure 6.1. The physical structure undergoes motion primarily in 
the Z direction, so all of the X and Y translations and Z rotations are con- 
strained in the finite element model. The two cross-beams at each end of the 
physical structure are connected by very short rods, and were originally mod- 
eled with NASTRAN ROD elements. However, this proved to be too stiff, so 
the rod elements were replaced with multi-point constraints constraining the 
Z-translational degree of freedom of nodes 9 and 10 to nodes 18 and 17, re- 
spectively. The payload simulator is supported by soft springs at nodes 11, 13, 
14, and 16 to simulate a bungee cord suspension system. Table 6.1 lists the 
undamped natural frequencies of the finite element model. The three lowest 
modes are the rigid body modes of the structure. 

The three interface degrees of freedom for this structure are the Z- 
translational degrees of freedom at nodes 4, 8, and 18. For the simulations 
presented in this section active excitation was used at all of the interface degrees 
of freedom, that is, no reaction forces were measured 

All of the simulations were performed in MATLAB© [65] running on 
a workstation. The mass and stiffness matrices were generated by and output 
from NASTRAN, and were then converted to a binary “MAT-file” readable by 
MATLAB©. To simulate viscous damping in the structure, modal damping 
at a level of 2% was added to all of the modes to obtain a physical damping 
matrix for the finite element model. 
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Tbble 6.1: Undamped Natural Frequencies of the Payload Simulator 


Number 

Frequency 

(Hz) 

Number 

Frequency 

(Hz) 

1* 

1.1564 

27 

1690.4 

2* 

6.0642 

28 

1702.2 

3* 

10.783 j 

29* 

2029.1 

4* 

18.127 

30* 

2031.2 1 

5* 

109.49 

31 1 

2079.9 

6* 

127.14 

32 

2104.3 

7* 

137.82 

33 

2272.3 

8* 

156.87 

34 

2289.9 

9* 

297.66 

35 

2746.6 ”j 

10* 

304.47 

36 

2750.7 

11 

512.54 

37 

2840.2 ] 

12 

521.24 

38 

2842.6 

13 

567.89 

39 

2943.9 

14 

570.46 

40 

2945.1 

15 

570.46 

41 

3403.3 

16 

583.43 

42 

3406.0 

17* 

606.45 

43 

5017.1 

18* 

631.08 

44 

5018.5 

19 

1073.0 

1 45 

7060.1 

20 

1078.6 

46 

7061.7 

21* 

1190.4 

1 47 

8887.0 

22* 

1193.3 

48 

8888.5 

23 

1341.5 

49 

11972. 

24 

1342.4 

50 

11973. 

25 

1342.4 

I 51 

14727. 

26 

1342.6 

52 

14727. 


^Target Modes — modes dominated by Z-translations 
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Figure 6.1: Payload Simulator Finite Element Model 


To generate the “experimental data” used in the simulations, the fre- 
quency response functions were generated by solving Eq. 3.6 for A(uk) given M, 
C, K, and F(uk). Except as noted, the input forcing function at each frequency 
u>k consisted of three independent unit forces applied in the Z-translational di- 
rection at nodes 4, 8, and 18 respectively. This data was then used as the input 
to the SSID algorithm to identify M, C, and K from the simulated measured 
acceleration responses, and force inputs, F(u>k). The input data used 

for the simulations presented in sections 6.2.4 and 6.3 contains additive noise 


while the input data for the other simulations is noise-free. 
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6.2 Simulation Results 

6.2.1 Identification of the Full-Order Model 

The first simulation demonstrates the algorithm s ability to identify 
structural models without any corrupting effects of noise, or of frequency trun- 
cation or spatial truncation. In order to identify all 52 modes, it was necessary 
to apply both forces and moments since the payload simulator is capable of de- 
forming in translational and rotational directions. The inputs were F z at node 
4, M x at node 17, and M„ at 18. The identification of the full-order model was 
based on 128 evenly spaced frequency lines ranging from 1 Hz to 15000 Hz. 
Computer memory constraints prohibited a finer frequency resolution. 

The results of the full-order identification are shown in Figs. 6.2-6.6. 
The individual entries of the estimated matrices were compared to those of 
the original matrices and the results are shown in Figs. 6.2-0.4. The SSID 
algorithm successfully identified the mass, damping, and stiffness matrices of 
the original structure to within 0.25%. The undamped natural frequencies were 
computed from the identified mass and stiffness matrices, and all 52 natural 
frequencies of the original structure were identified to within 8-digit accuracy. 
Figure 6.5 is a comparison between the exact drive point FRF at node 4 and 
one synthesized using the resulting identified M, C, and K matrices; Fig. 6.6 
shows the difference between the two FRFs. 

6.2.2 Identification of Reduced-Order Models 

In the previous simulation, the SSID algorithm was successfully able 
to identify the system matrices of the full-order 52 degree of freedom model 
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when all of the modes of the system were included. In a test environment, this 
would require response measurements at all 52 degrees of freedom, including 
the rotational degrees of freedom. However, it is impractical and impossible to 
measure the response of all of the degrees of freedom of a continuous structure. 
Therefore, spatial truncation is inevitable. 

The following simulations demonstrate the ability of the SSID algo- 
rithm to predict the dynamic characteristics of a structure when only a limited 
number of accelerometers, or other appropriate output devices, are available. 
The resulting identified model is referred to as a reduced-order model. Three 
different reduced-order models where used in the simulations. The number 
of degrees of freedom, or measurement locations, used for the reduced-order 
models was 10, 12 and 16 degrees of freedom. 

It is desired that the reduced-order models be able to represent the 
target modes of interest. The target modes to be identified are the rigid-body 
modes and other “global” modes of the structure (i.e., modes where localized 
bending predominates or where rotational DOFs predominate are excluded). 
These modes are indicated by an asterisk in Table 6.1. Note that the target 
modes are not just the lowest consecutive modes of the system, but are modes 
which are spaced throughout the spectrum. 

The Effective Interface method [66] was used for the process of se- 
lecting which degrees of freedom were to be “measured,” or retained, during 
the numerical simulations. The Effective Interface Method is an iterative sen- 
sor placement process that maximizes the observability of the target modes 
by maximizing the determinant of the Fisher information matrix. The Fisher 
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information matrix is formed by 

F = $l a *ta ( 6 - 1 ) 

where are the target modes partitioned to the retained degrees of freedom. 
The method examines the effect that the removal of each degree of freedom 
would have on the Fisher information matrix and removes one degree of freedom 
at a time. Table 6.2 summarizes the nodes selected by the Effective Interface 
Method for use as response locations for the three different reduced-order model 
sizes. In each case, the Effective Interface Method selected the Z-translational 
degree of freedom at the nodes listed in Table 6.2. 


Table 6.2: Reduced-Order Model Measurement Sensor Locations 


Model Size 

Measurement Locations 

10 

1, 2, 3, 4, 5, 6, 7, 8, 17, 18 I 

12 

1, 4, 5, 8, 11, 12, 13, 14, 15, 16, 17, 18 

16 

1, 2, 3, 4, 5, 6, 7, 8, 11, 12, 13, 14, 15, 16, 17, 18 


6.2.3 Effect of Input Frequency Spectrum on the SSID Algorithm 

For each of the reduced-order models, three different input frequency 
ranges were used to generate test data. This was done to determine how the 
SSID algorithm handles residual information. As was mentioned previously, 
the identified model depends on the frequency range. For each model size, 
better results can be obtained by adjusting the frequency range until the best 
frequency range is found. In a test environment, the “best” frequency range 
is unknown, so three input frequency ranges were arbitrarily selected. The 
different frequency ranges are summarized in Table 6.3. For each frequency 
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range, the frequency lines are equally spaced between the minimum frequency 
and the maximum frequency. The third input frequency range uses 1024 fre- 
quency lines instead of 512 since the frequency spacing using 512 would have 
been approximately 4 Hz and was considered too large to identify modes space 
2-3 Hz apart, which modes 21-22 and 29-30 are. 


Table 6.3: Input Frequency Spectrums 


Number 

Minimum 
Frequency (Hz) 

Maximum 
Frequency (Hz) 

Number of 
Frequency Lines 

HI 

1 

1 

300 

512 

0.584 

2 

1 

650 

512 

1.267 

3 

1 

2100 

1024 

2.050 


Table 6.4 compares the estimated undamped natural frequencies of 
the 10-DOF model for the different input frequency ranges. The estimated 
damping factors are given in Table 6.5. The damping factors are computed 
from the solution of the state space eigenproblem of Eq. 3.12. For the 10-DOF 
model, the best results are obtained for the frequency range of 1-300 Hz, which 
encompasses the first nine modes of the original system. As the maximum 
frequency is increased to include more modes, the resulting models tend to 
overestimate the natural frequencies and do a much poorer job estimating the 
damping factors, especially the lowest one. 

The resulting drive-point FRFs at node 4 from each 10-DOF model 
are shown in Figure 6.7. For the input frequency range of 1-300 Hz, the 
estimated response is almost identical to that of the response of the original 52 
degree of freedom model. The response of the identified model corresponding 
to the other two input frequency ranges overestimates the response near the 
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first mode. This is presumably caused by response of the higher modes, which 
the 10-DOF model is unable to represent. The 10-DOF model is trying to 
“fit” the response of the higher modes and thus overestimates near the first 
mode. Additionally, these two models do not estimate the amplitude of the 
true response near the resonances and anti-resonances as well as the first case 
does. 

In the FRF comparisons included in this report, the damping matrix 
used to form the synthesized FRFs was identified from the appropriate parti- 
tion of A a in Equation 3.15. Another damping matrix could also be formed by 
using the estimated damping factors and the real normal modes from the iden- 
tified mass and stiffness matrices. For a given identified model, this “modal” 
damping matrix generally overestimated the response at high frequencies while 
the identified damping matrix from Eq. 3.15 underestimates the response. A 
representative example is shown in Figure 6.8. 

To determine how the identified eigenvectors relate to the original 
ones, the Modal Assurance Criteria (MAC) [5] was computed. The MAC is 
not a true orthogonality check but is a scalar quantity that indicates the degree 
of independence between two mode shapes. The first 30 exact eigenvectors were 
partitioned to the retained degrees of freedom and compared to the eigenvectors 
obtained from the eigensolution of the identified mass and stiffness matrices. 
For all three 10-DOF models, the identified eigenvectors had MAC values of 
0.99 or greater for the corresponding target mode and values of less than 0.05 
for the remaining eigenvectors. 

The estimated undamped natural frequencies and damping factors 
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Table 6.4: Estimated Natural Frequencies — 10-DOF Model vs Input Fre- 
quency Spectrum 


Exact 

Fmax=300 Hz 

Fmax=650 Hz 

Fmax 

=2100 Hz 

Freq. 

Freq. 

Percent 

Freq. 

Percent 

Freq. 

Percent 

(Hz) 

(Hz) 

Error 

(Hz) 

Error 

(Hz) 

Error 

1.1564 

1.1552 

1.116e-01 

1.1301 

2.281e+00 

1.1453 

9.633e-01 

606.42 

6.0647 

6.755e-03 

6.0648 

9.276e-03 

6.0786 

2.359e-01 

10.783 

10.784 

7.202e-03 

10.956 

1.604e+00 

10.924 

1.305e+00 

18.127 

18.126 

5.608e-03 

18.133 

3.268e-02 

18.462 

1.850e+00 

109.49 

109.49 

3.624e-03 

109.60 

1.032e-01 

109.44 

4.352e-02 

127.14 

127.12 

1.203e-02 

127.02 

9.223e-02 

128.04 

7.052e-01 

137.82 

137.83 

4.228e-03 

137.84 

1.609e-02 

137.73 

6.875e-02 

156.87 

156.87 

1.263e-03 

156.78 

5.782e-02 

156.88 

7.647e-03 

297.66 

297.52 

4.736e-02 

296.70 

3.230e-01 

308.08 

3.502e+00 

304.47 

303.74 

2.376e-01 

304.12 

1.127e-01 

309.89 

1.781e+00 


Table 6.5: Estimated Damping Factors — 10-DOF Model vs Input Frequency 
Spectrum 


Exact 

Damp. 

Ratio 

Fmax=300 Hz 

Fmax=650 Hz 

Fmax=2100 Hz 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

0.02 

0.01974 

1.291e+00 

0.02674 

3.369e+01 

0.04127 

1.063e+02 

0.02 

0.01987 

6.368e-01 

0.01902 

4.921e+00 

0.09362 

3.681e+02 

0.02 

0.01895 

5.240e+00 

0.01163 

4.183e+01 

0.06398 


0.02 

0.01971 

1.469e+00 

0.01919 

4.028e+00 

0.07168 

2.584e+02 

0.02 

0.02000 

1.109e-02 

0.02010 

5.074e-01 

0.02017 

8.340e-01 

0.02 

0.01992 

3.767e-01 

0.01951 

2.458e+00 

0.02673 

WBMEBL 

0.02 

0.02000 

1.677e-02 

0.01992 

3.855e-01 

0.02126 

6.288e+00 

0.02 

0.01997 

1.290e-01 

0.02034 

1.724e+00 

0.02180 

9.002e+00 

0.02 

0.01994 

2.949e-01 

0.02227 

1.133e+01 

0.00683 

6.587e+01 

0.02 

0.01970 

1.474e+00 

0.02033 

1.644e+00 

0.00714 

6.432e+01 
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of the 12-DOF model for the different input frequency ranges are listed in 
Tables 6.6 and 6.7. The corresponding FRFs are shown in Figure 6.9. For the 
12-DOF model, the best results are obtained for the frequency ranges 1-300 Hz 
and 1-650 Hz. Note that the model corresponding to the input frequency range 
of 1-300 Hz was able to correctly estimate the system’s response for modes 11 
and 12 even though they are above the 300 Hz upper limit of the FRF data. 
The synthesized FRFs for these two models match very closely to those of the 
exact response, except near the response of the 12th mode. Here, the model 
corresponding to the first frequency range underestimates the response while 
the model corresponding to the 2nd input frequency range overestimates the 
magnitude of the response. 

The identified model corresponding to the third input frequency range, 
which includes all 16 target modes, overestimates the 11th and 12th natural 
frequencies and also does a much poorer job of estimating the damping factors. 
Again, the response of this model overestimates the response around the first 
mode and does a much poorer job of estimating response near the resonances 
and anti-resonances than do the models corresponding to the first two input 
frequency ranges. MAC values were also computed for the 12-DOF models 
with the same results as before, that is, near-perfect correlation between the 
identified modes and the target FEM modes. 

The estimated undamped natural frequencies and damping factors of 
the 16- DOF model for the different input frequency spectrums are listed in 
Tables 6.8 and 6.9, and the corresponding FRFs are shown in Figure 6.10. For 
the first 10 natural frequencies and damping factors, the identified model oorre- 
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Table 6.6: Estimated Natural Frequencies — 12-DOF Model vs Input Fre- 
quency Spectrum 


Exact 

Fmax=300 Hz 

Fmax=650 Hz 

Fmax 

=2100 Hz 

Freq. 

(Hz) 

Freq. 

(Hz) 

Percent 

Error 

Freq. 

(Hz) 

Percent 

Error 

Freq. 

(Hz) 

Percent 

Error 

1.1564 

1.1564 

1.465e-03 

1.1563 

1.560e-02 

1.1577 

1.073e-01 

6.0642 

6.0642 

3.454e-05 

6.0643 

6.787e-04 

6.0620 

3.767e-02 

10.783 

10.783 

8.332e-06 

10.783 

2.530e-03 

10.785 

1.427e-02 

18.127 

18.127 

1.780e-06 

18.127 

2.632e-03 

18.112 

7.970e-02 

109.49 

109.49 

5.694e-04 

109.49 

2.860e-03 

109.54 

4.366e-02 

127.14 

127.14 

1.377e-04 

127.19 

4.276e-02 

127.30 

1.261e-01 

137.82 

137.82 

7.153e-05 

137.82 


137.83 

5.248e-03 

156.87 

156.86 

8.364e-05 

156.86 

7.596e-03 

156.87 

4.080e-03 

297.66 

296.64 

6.828e-03 

296.99 

■MEMgsgllJ 

299.29 

5.461e-01 

304.47 

304.58 

3.859e-02 

304.41 

1.833e-02 

305.56 

3.598e-01 

606.45 

602.08 

7.211e-01 

605.93 

8.683e-02 

616.79 

1.704e+00 

631.08 

614.86 

2.570e+00 

609.92 

3.353e+00 

643.57 

1.979e+00 


sponding to the 1-300 Hz input frequency range performed the best. However, 
this model does not estimate the last four natural frequencies and damping val- 
ues very well. As shown by the MAC values in Table 6.10, the last two modes 
do not correspond to the any of the target modes. The poor estimation of the 
higher frequencies is also indicated in the corresponding synthesized FRF in 
Figure 6.10. 

The 16-DOF model corresponding to the second input frequency spec- 
trum estimates the first 12 frequencies very well. However, this model also 
identifies a frequency at 582.31 Hz, which agrees well with the 15th mode of 
the original 52 degrees of freedom model. This mode shape is similar to the 
2nd anti-symmetric bending mode of the structure, mode 9, except that the 
maximum displacement in the mode shape occurs in the X-rotational degree of 
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Table 6.7: Estimated Damping Factors — 12-DOF Model vs Input Frequency 
Spectrum 


Exact 

Damp. 

Ratio 

Fmax=300 Hz 

Fmax=650 Hz 

Fmax=2100 Hz 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

0.02 

0.02000 

9.909e-03 

0.02006 

2.965e-01 

0.01782 

1.088e+01 

0.02 

0.02000 

2.023e-03 

0.01995 

2.550e-01 

0.01854 

7.309e+00 

0.02 

0.02000 

4.948e-03 

0.02025 

1.257e+00 

0.01558 

2.209e+01 

0.02 

0.02000 

9.197e-04 

0.01936 

3.193e+00 

0.01929 

3.524e+00 

0.02 

0.02000 

1.234e-02 

0.02003 

1.563e-01 

0.01985 

7.721e-01 

0.02 

0.02000 

6.847e-03 

0.02012 

6.068e-01 

0.01855 

7.255e+00 

0.02 

0.02000 

1.396e-03 

0.01999 

4.555e-02 

0.01982 

8.789e-01 

0.02 

0.02000 

5.182e-03 

0.02000 

4.095e-03 

0.02000 

5.835e-03 

0.02 

0.01999 

7.062e-02 

0.01982 

8.829e-01 

0.03171 

5.855e+01 

0.02 

0.02004 

2.125e-01 

0.01989 

5.649e-01 

0.02753 

3.764e+01 

0.02 

0.02003 

1.565e-01 

0.01979 

1.029e+00 

0.018071 

9.641e+00 

0.02 

0.02003 

1.735e-01 

0.02774 

3.869e+01 

0.02286 

1.432e+01 


freedom at nodes 9, 10, 17, and 18. The inclusion of this mode is the suspected 
reason for the error in the estimation of the last 4 natural frequencies and for 
the error in the synthesized FRF for frequencies greater than 700 Hz. 

The identified model corresponding to the input frequency range of 
1-2100 Hz represents the dynamic characteristics of the original structure the 
best overall. All 16 natural frequencies are identified to within 1.0% and all of 
the MAC values are near unity for the corresponding target modes and less than 
0.05 for the remaining ones. The synthesized FRF matches the exact response 
almost exactly, except at very high frequencies, where it underestimates the 
magnitude of the true response. 

As indicated by the preceding simulations, the SSID-identified ma- 
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trices axe frequency-dependent. The best estimation of a given model order 
is when the input frequency range contains only the corresponding number of 
modes. For this reason, in the remaining simulations, the 10-DOF model was 
estimated using the first input frequency spectrum, the second input frequency 
range was used to estimate the 12-DOF model, and the 1-2100 Hz input fre- 
quency spectrum was used to estimate the 16-DOF model. 


Table 6.8: Estimated Natural Frequencies — 16-DOF Model vs Input Fre- 
quency Spectrum 


Exact 

Fmax=300 Hz 

Fmax=650 Hz 

Fmax 

=2100 Hz 

Freq. 

Freq. 

Percent 

Freq. 

Percent 

Freq. 

Percent 

(Hz) 

(Hz) 

Error 

(Hz) 

Error 

(Hz) 

Error 

1.1564 

1.1564 

5.902e-05 

1.1564 

6.472e-03 

1.1507 

4.948e-01 

6.0642 

6.0642 

3.468e-06 

6.0642 

4.260e-06 

6.0628 

2.430e-02 

10.783 

10.783 

3.174e-08 

10.794 

1.013e-01 

10.807 

2.165e-01 

18.127 

18.127 

3.659e-08 

18.124 

1.280e-02 

18.115 

6.439e-02 

109.49 

109.49 

1.992e-05 

109.48 

5.737e-03 

109.49 

1.539e-03 

127.14 

127.14 

1.238e-07 

127.16 

1.745e-02 

127.15 

7.669e-03 

137.82 

137.82 

2.793e-05 

137.82 

2.366e-04 

137.82 

2.866e-03 

156.87 

156.87 

4.226e-05 

156.86 

4.682e-03 

156.86 

5.341e-03 

297.66 

297.66 

1.745e-06 

296.70 

3.218e-01 

297.61 

1.796e-02 

304.47 

304.47 

2.041e-05 

304.47 

7.440e-04 

304.41 

1.745e-02 

606.45 

606.15 

5.092e-02 

582.31 

— 

606.24 

3.595e-02 

631.08 

613.54 

2.779e+00 

606.53 

1.319e-02 

626.89 

6.648e-01 

1190.4 

1047.1 

1.204e+01 

625.76 

8.430e-01 

1189.4 

8.538e-02 

1193.3 

1159.3 

2.852e+00 

1023.9 

1.399e+01 

1192.3 

8.922e-02 

2029.1 

1526.7 

2.476e+01 

1154.6 

3.243e+00 

2017.0 

5.990e-01 

2031.2 

1682.3 

1.718e+01 

1773.3 

1.261e+01 

2017.6 

6.686e-01 


For each of the reduced-order models, a “pseudo” drive-point FRF 
was estimated and compared to the original model. Using the identified system 
matrices, a drive-point FRF was estimated even though no force was actually 
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Table 6.9: Estimated Damping Factors — 16-DOF Model vs Input Frequency 
Spectrum 


Exact 

Fmax= 

=300 Hz 

Fmax= 

=650 Hz 

Fmax= 

=2100 Hz 

Damp. 

Ratio 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

0.02 

0.02000 

4.237e-04 

0.01997 

1.263e-01 

0.01980 

1.022e+00 

0.02 

0.02000 

4.263e-05 

0.02000 

7.409e-03 

0.02104 

5.179e+00 

0.02 

0.02000 

3.668e-06 

0.01859 

7.034e+00 

0.01824 

8.789e+00 

0.02 

0.02000 

6.994e-06 

0.01893 

5.333e+00 

0.02055 

2.754e+00 

0.02 

0.02000 

7.646e-04 

0.01999 

6.403e-02 

0.02009 

4.696e-01 

0.02 

0.02000 

3.306e-06 

0.01884 

5.778e+00 

0.01882 

5.911e+00 

0.02 

0.02000 

8.368e-05 

0.02000 

4.057e-03 

0.01993 

3.678e-01 

0.02 

0.02000 

5.035e-04 

0.01999 

6.897e-02 

0.01993 

3.481e-01 

0.02 

0.02000 

1.471e-04 

0.01966 

1.682e+00l 

0.02157 

7.847e+00 

0.02 

0.02000 H 

2.141e-04 

0.02000 

1.538e-03 

0.01980 

9.997e-01 

0.02 

0.02000 

2.160e-02 

0.02089 

4.473e+00 

0.02018 

8.870e-01 

0.02 

0.01996 

2.145e-01 

0.02002 

1.099e-01 

0.02593 

2.964e+01 

0.02 

0.01775 

1.126e+01 

0.02698 

3.492e+01 

0.0203 lT 

1.552e+00 

0.02 

0.02034 

1.185e+00 

0.02041 

2.049e+00 

0.01923 

3.830e+00 

0.02 

0.02135 

6.742e+00 

0.02017 

8.743e-01 

0.01549 

2.257e+01 

0.02 

0.02094 

4.685e+00 

0.02716 

3.580e+01 

0.01571 

2.145e+01 
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Table 6.10: Maximum MAC Values for 16-DOF Model vs Input Frequency 
Spectrum 


TAM 

Mode 

Fmax=300 Hz 

Fmax=650 Hz 

Fmax=2100 Hz 

FEM 

Mode 

Max 

MAC 

FEM 

Mode 

Max 

MAC 

FEM 

Mode 

Max 

MAC 

1 

1 

1.0000 

1 

1.0000 

1 

1.0000 

2 

2 

1.0000 

2 

1.0000 

2 

1.0000 

3 

3 

1.0000 

3 

1.0000 

3 

1.0000 

4 

4 

1.0000 

4 

1.0000 

4 

1.0000 

5 

5 

1.0000 

5 

1.0000 

5 

1.0000 

6 

6 

1.0000 

6 

1.0000 

6 

1.0000 

7 

7 

1.0000 

7 

1.0000 

7 

1.0000 

8 

8 

1.0000 

8 

1.0000 

8 

1.0000 

9 

9 

1.0000 

9 

1.0000 

9 

1.0000 

10 

10 

1.0000 

10 

1.0000 

10 

1.0000 

11 

18 

0.9987 

15 

0.8721 

17 : 

1.0000 

12 

17 

0.9344 

17 

1.0000 

18 

0.9970 

13 

20 

0.6239 

18 

0.9552 

21 

0.9999 

14 

21 

0.6639 

28 

0.6244 

22 

0.9999 

15 

32 

0.7196 

21 

0.8891 

30 

0.9999 

16 

27 | 

0.6063 

27 

0.3830 

29 

0.9992 
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applied at the degree of freedom during the identification of the model. The 
results of the comparisons are shown in Figures 6.11-6.13. As can been seen, 
good agreement is obtained between the estimated and the exact pseudo drive- 
point FRF. This suggests that interface information could possibly be obtained 
without applying a force or measuring a reaction force at the interface degrees 
of freedom. 

6.2.4 Effect of Noise on the SSID Algorithm 

In the simulations presented in the Section 6.2.1 through 6.2.3, the 
exact 52-DOF FRF data was used. The only “errors” that were introduced 
in the SSID processing of the exact FRF data was due to spatial and/or fre- 
quency truncation of the measured data. However, measured FRF data always 
contains some type of random error, or noise. The noise, for example, could 
be due to transducer error, signal processing and conditioning error, or other 
errors and uncertainties present in the measurement process. Therefore, for 
an identification algorithm to be applicable to test data, it must be robust 
enough to handle less than perfect data. The main concern of the simulations 
presented in this section is the numerical stability of the SSID algorithm in the 
presence of noise. 

The noise added to the simulated FRF data was uniformly distributed, 
pure random, and zero mean. A noise level of 2% was used in the simulations in 
this study. The noise level is the percent of the root-mean-square (RMS) value 
of the random noise to the RMS value of the noise free-signal. Acceleration is 
usually measured as the output of the structure, so the magnitude of the noise 
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is proportional to the magnitude of each accelerance FRF spectrum, not the 
receptance frequency response. A random phase error was also introduced; the 
maximum error on the phase was 2°. A typical “measured” response is shown 
in Figure 6.14. To reduce the effect of the noise, signal averaging was used, 
just as averaging would be employed in an actual test. The measured FRFs 
were averaged over 40 samples. 

The three least squares methods discussed in the previous Chapter 
were used to solve the over-determined system of equations, Eq. 3.7. No accept- 
able solution for any of the reduced-order models was ever obtained using the 
ordinary least squares method. The standard MATLAB least squares solver, 
which performs a Q-R decomposition, was used to solve the over-determined 
system of equations. The resulting identified system matrices were not positive 
definite, and typically yielded complex estimates for the undamped natural fre- 
quencies. Increasing the number of averages did not significantly improve the 
resulting solution. Figure 6.15 is a representative least squares solution and 
estimation of the drive-point FRF. This result indicates that in the implemen- 
tation of the SSID algorithm for use in test environments, a TLS solver will be 
necessary. 

A comparison of the estimations resulting from the different methods 
of solution for the 10-DOF model are shown Tables 6.11 and 6.12. Except for 
the first and the last mode, both methods of solution predict the undamped 
natural frequencies to within 1%. The STLS method gives slightly more accu- 
rate results, which is to be expected. For the 1st mode, which is where the noise 
affects the signal the most, the STLS estimation is much more accurate than 
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the TLS estimation. However, the STLS method does not accurately estimate 
the damping factors for the last three modes. 

The corresponding FRFs are shown in Figure 6.16. Again, both meth- 
ods of solution predict the response of the structure very well. The STLS 
method performs better than the TLS method in the low frequency response, 
while the TLS solution is more accurate than the STLS solution at higher fre- 
quencies. The TLS solution also overestimates the response near resonance for 
modes 4 and 5. 


Table 6.11: Estimated Natural Frequencies — 10-DOF Model vs Solution 
Method 


Exact 

Freq. 

(Hz) 

Noise Free 

TLS 

STLS 

Freq. 

(Hz) 

Percent 

Error 

Freq. 

(Hz) 

Percent 

Error 

Freq. 

(Hz) 

Peroent 

Error 

1.1564 

1.1552 

1.116e-01 

1.3706 

1.852e+01 

1.1795 

1.989e+00 

6.0642 

6.0647 

6.755e-03 

6.0654 

1.969e-02 

6.0633 

1.639e-02 

10.783 

10.784 

7.202e-03 

10.799 

1.434e-01 

10.780 

2.786e-02 

18.127 

18.126 

5.608e-03 

18.139 

6.756e-02 

18.118 

4.946e-02 

109.49 

109.49 

3.624e-03 

109.54 

4.068e-02 

109.55 

5.016e-02 

127.14 

127.12 

1.203e-02 

127.01 

1.040e-01 

127.05 

7.012e-02 

137.82 

137.83 

4.228e-03 

137.85 

2.018e-02 

137.85 

1.861e-02 

156.87 

156.87 

1.263e-03 

157.19 

2.003e-01 

157.33 

2.947e-01 

mm 

297.52 

4.736e-02 

299.96 

7.711e-01 

299.45 

6.012e-01 

304.47 

303.74 

2.376e-01 

312.42 

2.612e+00 

315.90 

3.754e+00 


Comparisons of the results for the 12-DOF model are shown Ta- 
bles 6.13 and 6.14 and in Figure 6.17. The same general trend is observed 
in this case as for the 10-DOF model. The STLS method provides a slightly 
better estimation of the dynamic characteristics of the original structure than 
the TLS solution for the low-frequency response. Both solutions overestimate 
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Table 6.12: Estimated Damping Factors — 10-DOF Model vs Solution Method 


Exact 

Damp. 

Ratio 

Noise Free 

TLS 

STLS 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

0.02 

0.01974 

1.291e+00 

0.01700 

1.502e+01 

0.02018 

9.057e-01 

0.02 

0.01987 

6.368e-01 

0.01990 

4.851e-01 

0.01994 

3.019e-01 

0.02 

0.01895 

5.240e+00 

0.01895 

5.226e+00 

0.01918 

4.099e+00 

0.02 

0.01971 

1.469e+00 

0.01967 

1.645e+00 

0.01978 

1.096e+00 

0.02 

0.02000 

1.109e-02 

0.01993 

3.698e-01 

0.01899 

5.025e+00 

0.02 

0.01992 

3.767e-01 

0.01876 

6.221e+00 

0.01921 

3.956e+00 

0.02 

0.02000 

1.677e-02 

0.02057 

2.860e+00 

0.01988 

5.937e-01 

0.02 

0.01997 

1.290e-01 

0.02319 

1.595e+01 

0.02265 

1.326e+01 

0.02 

0.01994 


0.01963 

1.841e+00 

0.01451 

2.744e+01 

0.02 

0.01970 

1.474e+00 

0.03474 


0.03621 

8.103e+01 


the response of the structure for frequencies greater than the input frequency 
range of 1-650 Hz. 

For the case of the 16-DOF model in the presence of noise, none of 
the three solution methods was able to provide an acceptable solution. The 
identified system matrices were not positive definite for any of the three meth- 
ods of solution. The erroneous solutions are more than likely due to the power 
to which the additive noise terms are raised in the frequency-domain equations 
of motion. This unsuccessful identification provides support to the claim that 
frequency-domain system identification methods work best as narrow-band es- 
timation procedures. Further evidence is provided by the successful identifica- 
tion employing the band processing approach, the results of which are listed in 


Table 6.17. 
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Table 6.13: Estimated Natural Frequencies — 12-DOF Model vs Solution 
Method 


Exact 

Noise Free 

TLS 

STLS 

Freq. 

Freq. 

Percent 

Freq. 

Percent 

Freq. 

Percent 

(Hz) 

(Hz) 

Error 

(Hz) 

Error 

(Hz) 

Error 

1.1564 

1.1563 

1.560e-02 

1.2983 

1.226e+01 

1.1701 

1.182e+00 

6.0642 

6.0643 

6.787en04 

6.0660 

2.812e-02 

605.85 

9.405e-02 

10.783 

10.783 

2.530e-03 

10.957 

1.605e+00 

10.830 

4.323e-01 

18.127 

18.127 

2.632e-03 

18.761 

3.500e+00 

18.221 

5.221e-01 

109.49 

109.49 

2.860e-03 

109.43 

5.631e-02 

109.46 

2.827e-02 

127.14 

127.19 

4.276e-02 

127.61 

3.710e-01 

127.49 

2.772e-01 

137.82 

137.82 

7.165e-04 

137.70 

8.719e-02 

137.73 

6.834e-02 

156.87 

156.86 

7.596e-03 

156.66 

1.359e-01 

156.71 

1.057e-01 

297.66 

296.99 

2.239e-01 

291.23 | 

2.160e+00 

295.51 

7.223e-01 

304.47 

304.4T 

1.833e-02 

304.49 

6.377e-03 

304.96 1 

1.610e-01 

606.45 

605.93 

8.683e-02 

600.29 

1.016e+00 

603.69 

4.558e-01 

631.08 

609.92 

3.353e+00 

653.34 

3.528e+00 

654. 151 

3.656e-KXT 


Table 6.14: Estimated Damping Factors — 12-DOF Model vs Solution Method 


Exact 

Noise Free 

TLS i 

5IL.£> 

1 

Damp. 

Ratio 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

r Damp. 
Ratio 

Percent 

Error 

0.02 

0.02006 

2.965e-01 

0.01737 

1.312e+01 

0.02002 

9.185e-02 

0.02 

0.01995 

2.550e-0i 

0.01975 

1.225e+00 

0.01989 

5.410e-01 

0.02 

0.02025 

1.257e+00 

0.02139 

6.977e+00 

0.02225 

1.127e+01 

0.02 

0.01936 

3.193e+00 

0.01758 

1.211e+01 

0.01717 

1.413e+01 

0.02 

0.02003 

1.563e-01 

0.02158 

7.932e+00 

0.02105 

5.268e+00 

0.02 

0.020i2 

6.068e-01 

0.01935 

3.234e+00 

0.01856 

7.178e+00 

0.02 

0.01999 

4.555e-02 

0.02171 

8.545e+00 

0.02154 

7.699e+00 

0.02 

0.02000^ 

4.095e-03 ' 

0.02102 

5.094e+00 

0.02089 

4.474e+00 

0.02 

0.01982 

8.829e-01 

0.02975 

4.876e+01 

0.02880 

4.400e+01 

0.02 

0.01989 

5.649e-01 

0.02251 i 

1.254e+01 

0.02257 

1.285e+01 

0.02 

0.01979 

1.029e+00 

0.02351 

1.754e+01 

0.01557 

2.215e+01 

0.02 

0.02774 

3.869e+01 

0.06257 

2.128e+02 

0.06822 

2.411e+02 
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6.3 Results of Narrow-Band Processing 

This section presents the results of simulations where the band pro- 
cessing approach described in Section 5.2 was used in the identification proce- 
dure. The band processing approach was motivated by the unsuccessful iden- 
tification of the 16-DOF model in the previous section. First, the method was 
applied to the 16-DOF model using noise-free data to determine if the idea 
would yield valid results. Then the data containing noise was used. 

The same simulated FRF data that was used in Section 6.2.3 for the 
16-DOF model, with the input frequency spectrum of 1-2100 Hz and 1024 
equally spaced frequency lines, is used here. Following the procedure outlined 
in Section 5.2, the frequency data was divided into 100-Hz frequency bands with 
50- Hz overlap, which results in 41 frequency bands. For each frequency band, 
the SSID identification proceeded through the solution of Eq. 3.12. The esti- 
mated natural frequencies computed using Eq. 5.19 are plotted in Figures 6.18- 
6.19. The vertical lines indicate the target-mode frequencies. Note that even 
though each frequency band spans only 100 Hz, the frequencies outside of a 
given band are also correctly identified. This is true only when noise-free data is 
used, as is evident when Figs. 6.18-6.19 and 6.21-6.22 are compared. However, 
the out-of-band frequencies will be discarded by virtue of the first selection 
criterion listed in Section 5.2. The roots in bands 1, 3, 6, 13, 24, and 40 were 
selected as the identified modal parameters, and the remaining portion of the 
SSID was then carried out to identify the system matrices. 

The three different approaches described in Section 5.2 for solving 
Eq. 3.23 were used to solve for the scale factors 5. These approaches are referred 
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to as BP1, BP2, and BP3 respectively. The BP1 estimates did not yield positive 
definite system matrices and axe not included here. The results from the other 
two methods are shown in Tables 6.15 and 6.16 and in Figure 6.20. The BP2 
es tima te provides the best estimation of the natural frequencies. However, 
the BP2 results tend to overestimate the response at high frequencies, while 
the response of BP3 is much closer to the one obtained using the original 
SSID formulation. The reason for the overestimation by the BP2 estimate 
has yet to be resolved. The estimates of the damping factors are the same 
for both band processing estimations, since the damping factors are computed 
from the same identified modal parameters. The difference between BP2 and 
BP3 is in the data used in Eq. 3.23 for the estimation of the fir coefficients. 
The band processing method estimates the damping factors very well. The 
results indicate that the modal parameters, Eq. 3.23, can be estimated band 
by band, collected and stored, and then used for the final estimation of the 
system matrices. 

The band processing approach was also applied to the noisy data 
for the 16 degree of freedom model used in the previous section. Using the 
original SSID formulation, no acceptable solution was ever obtained using any 
of the three least squares solution methods. The simulated data was divided 
into the same 41 frequency bands used in the last example and each band was 
processed through Eq. 3.12. The TLS method was used to solve Eq. 3.7. The 
estimated natural frequencies are plotted in Figures 6.21 and 6.22. Notice that 
when noise is applied to the signal, the TLS method tries to fit the data to the 
corresponding frequency band and most of the frequencies are estimated near 
the center of the band. This results in the identification of many “noise modes.” 
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Table 6.15: Estimated Natural Frequencies — 16-DOF Noise-Free Model Em- 
ploying Band Processing 


Exact 

Freq. 

(He) 

SSID 

BP2 

BP3 

Freq. 

(He) 

Percent 

Error 

Freq. 

(Hz) 

Percent 

Error 

Freq. 

(Hz) 

Percent 

Error 

1.1564 

1.1507 

4.948e-01 

1.1564 

2.972e-05 

1.1564 

1. 826e-05 

6.0642 

6.0628 

2.430e-02 

6.0642 

2.348e-04 

wm&M 

8.340e-04 

10.783 

10.807 

2.165e-01 

10.783 

1.096e-04 

10.783 

6.707e-06 

18.127 

18.115 

6.439e-02 

18.127 

9.541e-05 

18.126 

1.658e-03 

109.49 

109.49 

1.539e-03 

109.49 

3.223e-03 

109.49 

8.542e-04 


127.15 

7.669e-03 

127.14 

2.056e-03 

127.14 

1.502e-03 

137.82 

137.82 

2.866e-03 

137.81 

6.636e-03 

137.81 

9.261o-03 

156.87 

156.86 

5.341e-03 

156.87 

1.868e-03 

156.87 

4.088e-04 


297.61 

1.796e-02 

297.26 

1.345e-01 

297.56 

3.264e-02 

304.47 

304.41 

1.745e-02 

304.01 

1.495e-01 

304.50 

1.125e-02 

606.45 

606.24 

3.595e-02 

605.12 

2.202e-01 

606.46 

1.323e-03 

631.08 

626.89 

6.648e-01 

630.17 

1.450e-01 

629.75 

2.110e-01 

1190.4 


8.538e-01 

1188.9 

1.306e-01 

1195.1 

3.959e-01 


1192.3 

8.922e-02 

1191.1 

1.872e-01 

1196.4 

2.525e-01 

2029.1 

2017.0 

5.990e-01 

2029.0 

8.002e-03 

2144.5 

5.683e+00 

2031.2 

2017.6 

6.686e-01 

2031.8 

2.758e-02 

2171.2 

6.893e+00 
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Tbble 6.16: Estimated Damping Factors — 16-DOF Noise-Free Model Employ- 
ing Band Processing 


Exact 

Damp. 

Ratio 

SSID 

BP 

Damp. 

Ratio 

Percent 

Error 

Damp. 

Ratio 

Percent 

Error 

0.02 

0.01980 

1.022e+00 

0.02000 

2.688e-08 

0.02 

0.02104 

5.179e+00 

0.02000 

n&MMif 

0.02 

0.01824 

8.789e+00 

0.02000 

4.289e-09 

0.02 

0.02055 

2.754e+00 

0.02000 

1.311e-09 

0.02 


4.696e-01 

0.02000 

4.131e-06 

0.02 

0.01882 

5.911e+00 

E2E&3 

2.508e-08 

0.02 

0.01993 

■tgftisgjy 

0.02000 

7.951e-08 

0.02 

0.01993 

3.481e-01 

0.02000 

1.177e-07 

0.02 

0.02157 

7.847e+00 

0.02000 

3.463e-06 

0.02 

0.01980 

9.997e-01 

0.02000 

4.913e-08 

0.02 

EU 


0.02000 

1.115e-05 

0.02 

0.02593 

2.964e+01 

0.02000 

2.233e-02 

0.02 

0.02031 

1.552e+00 

0.02000 

1.140e-03 

0.02 

0.01923 

3.830e+00 

0.02000 

3.208e-05 

0.02 

0.01549 

2.257e+01 

0.02000 

1.718e-03 

0.02 

0.01571 

2.145e+01 

0.02000 

2.260e-03 
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An index of 0.95 was selected as the MPC criterion. Modes with MPC indices 
less than this value were rejected as noise modes. After application of the four 
selection criteria listed in Section 5.2, the selected roots are again plotted, as 
shown in Figures 6.23 and 6.24. 

The results of the identification are shown in Table 6.17 and in Fig- 
ure 6.25. The results are very encouraging, especially since no acceptable result 
was previously obtained using this same noisy data. As before, the BP2 solu- 
tion provides a better estimate of the natural frequencies. The high-frequency 
response of the BP2 estimation is an overestimation of the true response, while 
the low-frequency response of the BP3 estimate is an underestimation of the 
exact response. 


Table 6.17: Estimated Modal Parameters — 16-DOF Model Employing Band 
Processing with Noisy Data 


Exact 

Freq. 

(Hz) 

BP2 

BP3 

Exact 

BP 

Damp. 

Ratio 

Percent 

Error 

ESI 


ESI 


Damp. 

Ratio 

1.1564 

1.1435 

1.116e+00 

1.4273 

2.342e+01 

0.02 

0.00764 

6.178e+01 

6.0642 

6.0579 

1.054e-01 

6.1922 

2.110e+00 

0.02 

0.02016 

7.789e-01 

10.783 

10.784 

6.916e-03 

12.538 

1.629e+01 

0.02 

0.02046 

2.296e+00 

18.127 

18.132 

2.896e-02 

19.005 

4.846e+00 

0.02 

0.02143 

7.152e+00 

109.49 

109.51 

1.700e-02 

109.51 

1.989e-02 

0.02 

0.01880 

5.985e+00 

127.14 

127.24 

7.661e-02 

127.39 

2.013e-01 

0.02 

0.01977 

1.168e+00 

137.82 

137.82 

4.082e-04 

137.80 

1.642e-02 

0.02 

H 

4.474e-01 

156.87 

156.85 

1.089e-02 

156.87 

9.918e-04 

0.02 

0.01996 

1.836e-01 

297.66 

297.68 

5.516e-03 

297.92 

8.773e-02 

0.02 

0.01987 

6.248e-01 

304.47 

E3E1 

2.448e-02 

304.04 

1.415e-01 

0.02 

0.01975 

1.239e+00 

606.45 

usaa 

4.421e-02 


4.916e-02 

0.02 

0.02058 

2.911e+00 

631.08 

630.93 

2.382e-02 

631.76 

1.073e-01 

0.02 

0.01957 

2.145e+00 

1190.4 

1191.1 

5.449e-02 

1186.1 

3.625e-01 

0.02 

0.02004 

1.824e-01 

1193.3 

1193.1 

2.350e-02 

1192.0 

1.118e-01 

0.02 

0.02000 

1.637e-02 

2029.1 

2015.9 

6.522e-01 

2015.5 

6.693e-01 

0.02 

0.01864 

6.796e+00 

2031.2 

2031.8 

2.709e-02 

2054.4 

1.141e+00 

0.02 

0.01896 

5.213e+00 
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The band processing technique was also applied to the remaining two 
reduced-order models using noisy data. Figures 6.26 and 6.27 show the com- 
parisons of a drive-point FRF for each of the reduced-order models. The data 
used for the 12-DOF model was divided into 12 overlapping 100-Hz frequency 
bands. The results are similar to those obtained using the 16-DOF model. In 
order for the BP2 estimation to yield acceptable results, the 1 Hz data point 
was omitted. This would be similar to ignoring data with very low coherence 
values. The low frequency data is still severely corrupted by noise, even after 
sufficient averaging, and including this point in the BP2 estimation introduced 
large errors into the estimated frequency response. 

For the 10-DOF model, the BP2 estimation of the response is in seri- 
ous error, as indicated in Figure 6.27. For this model, the simulated data was 
split into 50- Hz frequency bands with 25- Hz overlap, resulting in 11 frequency 
bands. The reason for the error in the BP2 estimation is unknown. However, 
it is possible that the model is too small for use with the band processing 
technique or the frequency bands may too narrow. More research is needed 
to determine the width and the number of frequency bands to use for a given 
model order. The use or development of additional modal quality indicators 
should also be addressed to help identify and confirm the estimated structural 
modes. Also, coherence blanking or other measures might be useful in select- 
ing FRF data, or FRFs might be obtained by employing stepped sine-dwell 
processing. 
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6.4 Pseudo Degrees of Freedom 

This section presents the results of augmenting the model order using 
pseudo degrees of freedom (PDF). For this simulation, the 10-DOF model was 
used and the the input frequency range was 1-2100 Hz with 1024 equally space 
frequency lines. Using the pseudo degrees of freedom, the model order was 
expanded to create a 20-DOF pseudo-degree-of-freedom model, and the iden- 
tification process through Eq. 3.12 was completed. As mentioned in Section 
5.3, the remaining portion of the identification has yet to be defined for pseudo 
degrees of freedom. No noise was added to the data. Table 6.18 gives a listing 
of the estimated frequencies from the solution of Eq. 3.12 and compares those 
to the frequencies obtained using only the original ten degrees of freedom. 


Table 6.18: Identified Natural Frequencies Using Pseudo Degrees of Freedom 


Pseudo DOF Model 

] Original Model 

Mode 

R-eq. 

Mode 

BYeq. 

Mode 

Freq. 

Number 

(Hz) 

Number 

(Hz) 

Number 

(Hz) 

1 

1.1026 

11 

304.24 

1 

1.1454 

2 

6.0572 

12 

463.91 

2 

6.0878 

3 

10.387 

13 

497.10 

3 

10.989 

4 

18.063 

14 

1189.8 

4 

18.293 

5 

109.50 

15 

1193.2 

5 

109.52 

6 

128.19 

16 

1675.1 

6 

127.75 

7 

137.82 

17 

1686.4 

7 

137.81 

8 

156.89 

18 

1713.0 

8 

156.99 

9 

271.10 

19 

2029.5 

9 

308.15 

10 

297.17 

20 

2031.6 

10 

311.36 


Using the original 10 degrees of freedom, only the lowest 10 modes 
could be identified even though more are present in the data. When the order 
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of the model is expanded by use of pseudo degrees of freedom, the model is 
able to identify 4 additional target modes, modes 21, 22, 29, and 30, plus 
some additional computational modes. Some sort of modal quality indicator 
is needed to distinguish between the physical and computational modes of the 
PDF model. This approach doubles the model order. In some instances, this 
may be too large and the model order may need to be decreased from twice 
the original size. A reduction step would then be required to determine the 
appropriate model order. The initial results are encouraging, since some of the 
higher frequencies can be identified using the pseudo degrees of freedom. 

The mode shapes from the PDF model are listed in Tables 6.19 
and 6.20. The mode shapes listed here are the normal modes of the PDF 
model. They are obtained by solving the algebraic eigenvalue problem using 
/C = which is obtained from the solution of Equation. 3.7. Note that 

for the lower frequency modes, the mode shape appears in the lower-frequency 
partition of the mode shape and the higher frequency mode shapes appear in 
the high-frequency partition, or in the pseudo degrees of freedom. The pseudo 
degrees of freedom could possibly be treated as generalized coordinates and 
then the SSID algorithm could proceed with the identification of the system 
matrices. 
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Table 6.19: Mode Shapes of the 10-DOF Model Employing Pseudo Degree of 
Freedom - Modes 1 through 10 
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Table 6.20: Mode Shapes of the 10-DOF Model Employing Pseudo Degree of 
Freedom - Modes 11 through 20 


Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

Mode 

ii 

12 

13 

14 

15 

16 

17 

18 


20 

- 1.000 

0.240 

mmizm 



1.000 

-0.095 

0.095 

mwm 

0.018 

0.650 

- 1.000 

- 1.000 


Bl 

0.695 

-0.032 

0.032 

BBSS 

0.047 


0.760 

0.367 

-0.023 

EMU 

BBS 


0.046 

E!2!£2fl 


0.938 




-0.086 

0.839 

Biliiil 

0.087 

BE21 

beei 

- 1.000 

-0.240 


-0.079 

0.054 

- 1.000 

BBS 

0.095 

-0.042 

BEB1 

0.650 

1.000 

- 1.000 

-0.042 

0.030 


-0.032 

0.032 

BMi 

BBS 

-0.672 

-0.760 

0.367 

MSHStEM 

0.020 

0.147 

BfiaJifl 

0.046 

BEES 

mssM 

0.938 

0.814 

Bll:K:tl 


0.086 


gllttirJ 

0.087 

BBS 

0.013 

g»MHi 


0.282 

0.018 

0.000 

0.000 

0.074 


0.019 

0.000 

0.542 

0.000 

0.317 

0.020 

0.000 

0.000 

E!2!£il 

0.006 

0.006 

0.000 

0.000 

0.062 

0.015 



0.028 

BilM:l 

0.158 

0.420 

0.439 

0.000 

0.080 

0.034 

-1.000 

- 1.000 


E!ZU£!l 

0.080 

0.995 

0.995 

0.000 

gtlt/rl 

B*HK*>i 

0.999 

1.000 

0.050 

0.063 

BEES 

1.000 

1.000 

lEzoa 

rnmm 

0.446 

0.375 

0.384 

gilUrj 

0.186 

BBS 

0.436 

0.428 

0.000 

-0.062 

0.015 

-0.371 

0.384 1 

BEES 

gjwrei 

0.158 

0.420 


0.000 

-0.080 

0.034 

- 1.000 

1.000 

0.054 


0.080 

0.995 

-0.995 

0.000 

0.077 

Btllttfrl 

0.999 

- 1.000 


0.063 

boot*! 

1.000 

- 1.000 

IBE&il 

w»K.’:»T 

0.446 

0.375 


0.115 

0.186 


0.436 

Bl JEW 

0.000 

0.000 

WillH 


0.000 

0.000 

1.000 

- 1.000 

bees 

0.000 

0.011 

0.000 

0.289 

0.100 

0.000 

0.000 


0.846 

BBS 

0.000 
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6.5 Comparison of S SID-Identified Models with Exist- 
ing TAMs 

In this section, the SSID-identified reduced-order models are com- 
pared to three of the TAM models discussed in Chapter 4: the Guyan model, 
the Improved Reduction System (IRS) model, and the Craig-Bampton (C-B) 
model in physical coordinates. The appropriate fixed-interface normal modes 
for use in the C-B models were selected using the technique described in Section 
4.4.2. For each TAM, the estimated natural frequencies and the prediction of 
the dynamic response are compared for each of the reduced-order model sizes, 
10, 12, and 16 degrees of freedom. Orthogonality and cross-orthogonality checks 
were performed using each TAM, and very good results were found for each 
TAM, with diagonal values near unity and off-diagonal terms less than 0.1. 

Tables 6.21 through 6.23 list the estimated undamped natural fre- 
quencies for each of the reduced-order model sizes. Except for the very lowest 
frequencies, the SSID reduced-order models generally provide a better estima- 
tion of the natural frequencies than does the Guyan reduction method. The 
IRS and C-B models also estimate the lower natural frequencies better than 
the do SSID models. However, for the middle and higher frequencies of each 
model size, the SSID-identified model performs as well, if not slightly better, 
than these other two TAMs. 

A comparison of the estimated FRFs for each of the model sizes is 
given in Figures 6.28-6.30. The physical damping matrix used for each of the 
TAMs was computed by using the real normal modes estimated by each model 
and the damping factors estimated by the SSID algorithm for the corresponding 
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model size and computing a modal damping matrix. As shown in the figures, 
all four methods represent the dynamic response of the original 52-DOF model 
very well. 


Table 6.21: Estimated Natural Frequencies — 10-DOF TAMs 


1 SSID 

Guyan 

■ 

[RS ] 


3-B 


Percent 


Percent 


Percent 

Freq. 

Percent 

nil 

Error 


Error 

mSm 

Error 

(Hz) 

Error 

1.1552 

1.116e-01 

1.1564 

1.100e-04 

1.1564 

3.834e-08 

1.1564 

2.352e-06 

6.0647 

6.755e-03 

6.0644 

2.576e-03 

6.0642 

1.064e-09 

6.0643 

1.660e-04 

10.784 

7.202e-03 

10.785 

1.004e-02 

10.783 

1.411e-09 

10.783 

3.144e-C4 

18.126 

5.608e-03 

18.132 

2.629e-02 

18.127 

2.559e-09 

18.127 

2.187e-03 

109.49 

3.624e-03 

109.56 

6.129e-02 

109.49 

1.043e-05 

109.52 

2.261e-02 

127.12 

1.203e-02 

127.46 

2.558e-01 

127.14 

3.221e-04 

127.32 

1.400e-01 

137.83 

4.228e-03^ 

138.08 

1.866e-01 

137.82 

3.570e-05 

137.83 

5.044e-03 

156.87 

1.263e-63 

157.25 

2.436e-01 

156.87 

9.836e-04 

157.06 

1.194e-01 

297.52 

4.736e-02 

309.&T 

4.087e+00 

297.82 

5.248e-02 

300.04 

7.981e-01 

303.74 

2.376e-01 

314.79 

3.392e+00 

304.51 

1.368e-02l 

306.30 

6.034e-01 


Table 6.22: Estimated Natural Frequencies — 12-DOF TAMs 


1 SSID 1 

Guyan 

IRS 


C-B 

Freq. 

Percent 

Freq. 

Percent 

Freq. 

Percent 

freq. 

Percent 

(Hz) 

Error 

(Hz) 

Error 

(Hz) 

Error 

(Hz) 

Error 

1.1563 

1.560e-02 

1.1564 

4.802e-06 

1.1564 

3.876e-08 

1.1564 

4.419e-07 

6.0643 

6.787e-04 

6.0643 

8.648e-05 

6.0642 

1.052e-09 

6.0642 

2.380e-05 

10.783 

2.530e-03 

10.783 

6.218e-04 

10.783 

7.343e-10 

10.783 

2.519e-04 

18.127 

2.632e-03 

18.127 

2.137e-03 

18.127 

1.610e-09 

18.127 

1.752e-03 

109.49 

2.860e-03 

109.55 

5.700e-02 

109.49 

5.814e-06 

109.50 

5.780e-03 

127.19 

4.276e-02 

127.40 

2.055e-01 

127.14 

1.957e-04 

127.28 

1.120e-01 

137.82 

7.165e-04 

137.84 

1.035e-02 

137.82 

9.817e-07 

137.83 

4.813e-03 

156.86 

7.596e-03 

157.01 

8.613e-02 

156.87 

3.124e-05 1 

156.90 

1.518e-02 

296.99 

2.239e-01 

301.40 

1.257e+00 

297.76 

3.218e-02 

299.54 

6.327e-01 

304.41 

1.833e-02 

306.93 

8.094e-01 

304.48 

3.440e-03 

304.69 

7.437e-02 

605.93 

8.683e-02 

608.57 

3.492e-O^T 

595.82 

1.753e+00 

607.42 

1.589e-01 

609.92 

3.353e+00 

611.20 1 

3.151e+00 

606.46 

3.901e+00 

617.48 

2.155e+00 


Percent differences were computed between the SSID mass and stiff- 
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Table 6.23: Estimated Natural Frequencies — 16-DOF TAMs 


SSID 

Guyan 

IRS 

C-B 

Freq. 

Percent 

Freq. 

Percent 

Freq. 

Percent 

Freq. 

Percent 

(Hz) 

Error 

(Hz) 

Error 

(Hz) 

Error 

(Hz) 

Error 

1.1507 

4.948e-01 

1.1564 

1.864e-06 

1.1564 

3.904e-08 

1.1564 

1.242e-07 

6.0628 

2.430e-02 

6.0643 

7.634e-05 

6.0642 

1.108e-09 

6.0642 

5.057e-06 

10.807 

2.165e-01 

10.783 

3.467e-01 

10.783 

7.736e-10 

10.783 

2.255e-04 

18.115 

6.439e-02 

18.127 

2.036e-03 

18.127 

1.557e-09 

18.127 

1.569e-03 

109.49 

1.539e-03 

109.50 

8.567e-03 

109.49 

5.435e-06 

109.49 

1.934e-03 

127.15 

7.669e-03 

127.27 

1.026e-01 

127.14 

1.943e-04 

127.27 

1.006e-01 

137.82 

2.866e-03 

137.83 

7.427e-03 

137.82 

6.314e-07 

137.83 

3.774e-03 

156.86 

5.341e-03 

156.91 

2.329e-02 

156.87 

2.836e-05 

156.88 

5.222e-03 

297.61 

1.796e-02 

299.12 

4.909e-01 

297.75 

3.140e-02 

299.38 

5.779e-01 

304.41 

1.745e-02 

304.53 

2.095e-02 

304.47 

6.180e-06 

304.51 

1.448e-02 

606.24 

3.595e-02 

604.31 

3.535e-01 

595.50 

1.807e+00 

606.61 

2.544e-02 

626.89 

6.648e-01 

607.01 

3.814e+00 

606.43 

3.906e+00 

617.20 

2.200e+00 

1189.4 

8.538e-02 

1212.4 

1.848e+00 

1133.1 

4.816e+00 

1189.5 

7.938e-02 

1192.3 

8.922e-02 

1214.1 

1.738e-f00 

1190.8 

2.093e-01 

1191.1 

1.911e-01 

2017.0 

5.990e-01 

2227.5 

9.774e+00 

1752.9 

1.361e+01 

2029.3 

7.551e-03 

2017.6 

6.686e-01 

2228.7 

9.721e+00 

2052.1 

1.026e+00 

2035.0 

1.865e-01 


ness matrices and those of the other TAMs. Figures 6.31 and 6.32 are repre- 
sentative results from the different comparisons. The only difference is that 
for the 10-DOF models, the percent differences are slightly smaller and for the 
16-DOF models, the differences are slightly larger. The main point is that the 
mass and stiffness matrices of the various TAM models are noticeably differ- 
ent. They do, however, yield essentially the same natural frequencies and result 
in the same estimated frequency response functions. Berman points out that 
multiple solutions of the same model can exist and have significantly different 
terms in the system matrices [67]. He shows that the stiffness matrix can be 
written in terms of the modes of the system and that if not all of the modes 
of the system are included, the estimated stiffness matrix will not be the same 
as the original stiffness matrix, but will still be representative of the system’s 




stiffness over a specified frequency range. This leads to the conclusion that, in 
test situations, it is impossible to identify a unique mass and stiffness matrices. 
However, the identified matrices still represent the dynamic characteristics of 
the original structure. 

Since the SSID algorithm identifies reduced-order mass and stiffness 
matrices of the physical system, all the information associated with a Craig- 
Bampton model is available, including information about the constraint modes. 
The SSID reduced-order models were put into the Craig-Bampton format and 
then compared to the corresponding Craig-Bampton analytical models gener- 
ated from the 52-DOF FEM. Tables 6.24 through 6.26 compare the eigenvalues 
of the fixed-interface normal modes of the Craig-Bampton models to those ob- 
tained from the SSID-identified models in Craig-Bampton format. The data 
associated with the constraint modes is compared by examining the appropriate 
partition of the stiffness matrices and is shown in Figures 6.33 through 6.35. 
As indicated by the results, the SSID-identified models represent the fixed- 
interface and constraint-mode information very well. 
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Table 6.24: Comparison of the Eigenvalues of the Fixed-Interface Normal 
Modes for the 10-DOF Model 



(rad/sec 2 ) 




2.4449e+04 


3.3205e+05 


3.5364e+05 


6.9749e+05 


2.641 le+06 


2.6885e+06 



6.9753e+03 


2.4535e+04 


3.3086e+05 


3.5389e+05 


6.9720e+05 


2.6397e+06 


2.6608e+06 


Percent 

Difference 


2.039e-01 


3.511e-01 


3.602e-01 


7.031e-02 


4.111e-02 


5.427e~02 


1.033e+00 



Table 6.25: Comparison of the Eigenvalues of the Fixed-Interface Normal 
Modes for the 12-DOF Model 
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Table 6.26: Comparison of the Eigenvalues of the Fixed-Interface Normal 
Modes for the 16-DOF Model 


C-B 

SSID 


Ar 

A r 

Percent 

(rad/sec 2 ) 

(rad/sec 2 ) 

Difference 

6.9895e+03 

7.0152e+03 

3.682e-01 

2.4449e+04 

2.5629e+04 

4.825e+00 

3.3205e+05 

3.3724e+05 

1.560e+00 

3.5364e+05 

3.5690e+05 

9.217e-01 

6.9749e+05 

6.9884e+05 

1.944e-01 

2.641 le+06 

2.6201e+06 

7.957e-01 

2.6885e+06 

2.6923e+06 

1.412e-01 

1.2978e+07 

1.2986e+07 

6.287e-02 

1.4456e+07 

1.3386e+07 

7.404e+00 

5.3240e+07 

5.3379e+07 

2.603e-01 

5.3430e+07 

5.4037e+07 

1.138e+01 

1.5448e+08 

1.4891e+08 

3.608e+00 

1.5476e+08 

1.4980e+08 

3.204e+00 


6.6 SSID Implementation with Reaction Forces Included 

This section presents simulations when interface reaction forces were 
included in data used in the identification of a full-order model of the substruc- 
ture. In this simulation, reaction forces where measured at all three interface 
degrees of freedom. At each interface degree of freedom, a 9-DOF test stand was 
attached to the substructure, resulting in a 79 degree of freedom substructure- 
test stand assembly. The resulting coupled model is statically determinate. 
Two independent excitation forces, located at nodes 11 and 16, where used in 
the identification of the substructure. 


The original 52-DOF model was first reduced to the 16 Z-translational 
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degrees of freedom via Guyan reduction, and the resulting model was con- 
sidered to be the “true” substructure. The same was done for the coupled 
79-DOF substructure-test stand assembly, resulting in a 25-DOF model. Ta- 
ble 6.27 lists the natural frequencies of the substructure-test stand assembly 
and of the substructure alone. The dynamic response of the coupled system 
is significantly different from that of the substructure alone. Note that lowest 
natural frequency of the coupled system is significantly higher than that of the 
substructure alone. Figure 6.36 compares the response at the same degree of 
freedom for the substructure-test stand assembly (B structure) to that of the 
substructure alone (A structure). Rayleigh damping was used instead of modal 
damping to insure that there was no unwanted damping coupling between the 
substructure and the test stand. 

A frequency spectrum of 1-1000 Hz with 512 equally spaced frequency 
fines was used for the identification. The results of the simulation are shown in 
Figures 6.37 through 6.41. The results show that, with noise-free data, the SSID 
algorithm was able to successfully identify the substructure from the coupled- 
system data. The system matrices were all identified to less than 0.03%, and 
the predicted response matches the original response very well. The estimated 
undamped natural frequencies and damping factors were all identified to within 
0.01%. Note that the highest four frequencies were successfully identified even 
though they are well above the range of the FRF data used in the identification. 

Simulations of the full-order model with the reaction forces included 
were also attempted using noisy data. However, to date none has been suc- 
cessful. One possible reason is the way the noise was applied and the way the 
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Table 6.27: Undamped Natural Frequencies of the Coupled System and of the 
Substructure Alone 


B Substructure j 

A Substructure 

Number 

Frequency 

(Hz) 

Number 

Frequency 

(Hz) 

1 

13.304 

1 

1.1564 

2 

24.775 

2 

6.0643 

3 

90.682 

3 

10.783 

4 

93.526 

4 

18.127 

5 

132.47 

5 

109.50 

6 

253.12 

6 

127.27 

7 

254.95 

7 

137.83 

8 

539.25 

8 

156.91 

9 1 

541.59 

9 

299.12 

10 

769.29 

10 

304.53 

11 

782.15 

11 

604.31 

12 

1223.4 

12 

607.01 

13 

1224.6 

13 

1212.4 

14 

1277.3 

14 

1214.1 

15 

1415.2 

15 

2227.5 

16 

1416.0 

16 

2228.7 

17 

1561.1 



18 

2230.1 



1 19 

2231.6 



20 

4080.2" 



21 

4080.3 



22 

! 4081.4 



23 

6118.4 



24 

6118.4 



25 

6118.8 
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reaction force was computed. Recall that the noise is applied such that it is 
proportional to the magnitude of the accelerance FRF. For the coupled sys- 
tem, the acceleration response at low frequencies is very small, while at high 
frequencies it is several orders of magnitude larger. The opposite is true for the 
displacement response, and the reaction force was computed as the difference 
between the displacement of the substructure-test stand interface DOFs mul- 
tiplied by the coupling reaction-spring constant. Even after averaging, the low 
frequency acceleration response and the high frequency displacement response, 
and thus the measured reaction force, are still quite corrupted by the noise and 
these responses are believed to be very important in these simulations where 
the reaction forces are included. This indicates that very accurate measure- 
ment of the reaction forces will be required for use in the SSID algorithm. An 
example of a “measured” accelerance FRF of the coupled system is shown in 
Figure 6.42. 

The measured reaction force FRFs due to the force at node 11 that 
were used in the identification of the full-order, noise-free simulation are pre- 
sented in Figure 6.43. Note the low response at high frequencies for the in- 
terface forces at nodes 4 and 18. To alleviate the problems relating to the 
reaction forces, several different ideas, such as varying the stiffness of the inter- 
face springs, were attempted. Simulations using mass-loaded interfaces were 
also attempted to increase the magnitude of the measured reaction force at 
higher frequencies, but none has proved successful. 

Additional simulations have been attempted to identify a reduced- 
order model with the reaction forces included. At this time, no successful 
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identification of the “payload” used in the present study has been made, and 
this is still an active area of research. Successful reduced-order-model results 
have been obtained with the reaction forces included for lumped-parameter 
models; these can be found in Reference [63]. 
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Figure 6.6: Difference between the Exact and Identified Drive-Point FRFs 




Figure 6.7: Comparison 
Spectrum 


ffl 



Frequency (Hz) 
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Figure 6.8: Comparison of FRFs Based on Two Different Identified Damping 
Matrices 
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Figure 6.9: Comparison of FRFs for the 12-DOF Model vs Input Frequency 
Spectrum 
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Figure 6.10: Comparison of FRFs for the 16-DOF Model vs Input Frequency 
Spectrum 
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Figure 6.11: Comparison of Pseudo Drive-Point FRFs for 10-DOF Model 
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Figure 6.14: Typical “Measured” FRF with Added Noise 
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Figure 6.15: A Least Squares Solution in the Presence of Noise 
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Figure 6.16: Comparison of FRFs of 10-DOF Model vs Solution Method 
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Figure 6.18: Estimated Natural Frequencies of 16-DOF Noise-Free Model Em- 
ploying Band Processing 
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Figure 6.19: Estimated Natural Frequencies of 16-DOF Noise-Free Model Em- 
ploying Band Processing (continued) 
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Figure 6.20: Comparison of FRFs of 16-DOF Noise- Free Model Using Band 
Processing 
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Figure 6.21: Estimated Natural Frequencies: Band Processing Using Noisy 
Data 
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Figure 6.23: Estimated Natural Frequencies After Mode Selection: Band Pro- 
cessing Using Noisy Data 
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Figure 6.24: Estimated Natural Frequencies After Mode Selection: Band Pro- 
cessing Using Noisy Data (continued) 
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Figure 6.25: Comparison FRFs of 16-DOF Model Employing Band Processing 
with Noisy Data 
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Figure 6.27: Comparison FRFs of 10 DOF Model Employing Band Processing 
with Noisy Data 
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Figure 6.31: Comparison of Mass Matrices for the 12-DOF TAMs — IRS vs 
SSID 



Figure 6.32: Comparison of Stiffness Matrices for the 12-DOF TAMs — IRS 
vs SSID 
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Figure 6.33: Comparison of the Constraint-Mode Partition of the C-B Stiffness 
Matrix for the 10-DOF Model 



Figure 6.34: Comparison of the Constraint-Mode Partition of the C-B Stiffness 
Matrix for the 12-DOF Model 
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Figure 6.36: Comparison of FRFs for the Coupled System and for the Sub- 
structure Alone 
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Figure 6.37: Percent Error in the SSID-Estimated Mass Matrix 



Figure 6.38: Percent Error in the SSID-Estimated Damping Matrix 
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Figure 6.39: Percent Error in the SSID-Estimated Stiffness Matrix 
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Figure 6.42: Noisy Accelerance FRF of the Coupled System 







Chapter 7 


Conclusions and Recommendations 


The purpose of this investigation was to examine the performance of 
a new structural identification algorithm, the Substructure System Identifica- 
tion algorithm, using numerical simulations for a moderate size finite element 
model. The algorithm successfully identified the substructure when active ex- 
citation was used at all the interface degrees of freedom and also for the case 
when reaction forces where measured at all the interface locations. The ef- 
fects of spatial and frequency truncation were studied and it was found that 
the algorithm is capable of identifying valid reduced-order structural models. 
Simulations were also performed using data containing noise and it was shown 
that the algorithm is robust enough to handle such data. 

A technique referred to as band processing was presented and shown 
to successfully identify structural models from noisy frequency data having a 
broad bandwidth. Thus, the SSID algorithm can be used as both a narrow-band 
and broad-band frequency-domain identification procedure. Pseudo degrees of 
freedom were examined as a way to expand the model size when there are more 
modes present in the data than there are output sensors. The model employing 
pseudo degrees of freedom was able to identify some of the additional modes 
when the original model could not. 

The identified system matrices from the SSID algorithm were also 
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compared to other test-analysis models and it was found that the SSID reduced- 
order models represent the dynamic characteristics of the structure as well as, if 
not better than, the other test-analysis models. In addition, it was shown that 
the SSID reduced-order models provide the information necessary to obtain 
the fixed-interface modal data associated with the Craig-Bampton substructure 
model and the data associated with the constraint modes as well. 

The following recommendations are made regarding the SSID algo- 
rithm and future work in this area: 

• Remove the dependence on the scaling of the fixed-interface normal modes 
and determine an absolute measure that can be used to select the fixed- 
interface normal modes that contribute most to the target modes for the 
Craig-Bampton reduced-order models. 

• Employ additional simulations to examine the identification of reduced- 
order models when the measured the reaction forces are included. Also, 
research is needed in the area of the measurement of reaction forces to 
determine acceptable reaction force levels. 

• Investigate additional modal quality indicators for use in band process- 
ing and with models employing pseudo degrees of freedom to distinguish 
between structural and computational modes. 

• Explore the band processing technique further, such as determining the 
optimal bandwidth and frequency resolution for different model sizes and 
modal densities. In addition, the SSID algorithm using the band process- 
ing technique should be implemented on a parallel-computing machine. 
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• Further study the area of proper model order determination for the SSID 
algorithm. The necessary steps needed to use pseudo degrees of freedom 
to increase the model order in the second stage of the SSID algorithm 
should be determined. Also, examine the case when there are much fewer 
modes than sensors. In this case, the resulting SSID model will not be 
a minimal-order model, but will be larger and the order of the model 
should somehow be reduced. 



Appendix A 

MATLAB© Source Code 


% Set Parameters: 

% ne = No. Excitation Frequencies 
% nx = No. of DOF of Substructure Model 
% n f = No. of Active Forces 
% nr = No. of Reactions 
% ns = No. of States 

% nb = No. of DOF of Substructure and Test Stand 
% nsample = No. of samples of excitation source 
% 

clear rand(‘seed’,0); 

% 

na=size(KA,l); 

% Compute the mass-normalized eigenvectors 
fphia,eiga]=massnorm(KA,MA) ; 
freqa=sqrt (eiga) / (2*pij ; 

% Compute the modal damping matrix 
jCA,zetaa]=damping(MA,KA,ones(na,l)*.02); 

% Forced DOF: 10, 22, and 50 

% Reaction DOF: NONE 

fdof=[10 22 50]; nf=3; nr=0; nd=nf+nr; 

% 

% “Preserve” these modes of the “A” Structure 
modes=[l:10, 17, 18, 21,22, 29,30]; %16 DOF Model 

%modes=[l:10,17,18]; %12 DOF Model 

%modes=[l:10j; %10 DOF Model 

% DOF to RETAIN from the Original Model 
rdof=eidv(phia(: , modes) , [1 :na] 1 .length (modes)) ; 

% DOF to OMITT from the Original Model 
odof=omitdof(na,rdof); 

% Retained Interior DOF of the Original Model and the Forced DOF 
% of the Reduced Model 

ja,ridof,frdof] =keptdof (na,odof ,fdof ( 1 :nf ) ,fdof (nf+ 1 :nd)) ; 

% Size of the substructure to identify 
nx=length(rdof); ns=2*nx; 

% 
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% Input Test Frequencies 

ne=512; deltaf=l; fmin=1.0; finax=300; 

Freq=linspace(fmin , fmax ,ne) ; W=Freq*2*pi; W2=W. "2; 

% 

% Noise parameters 

nsample=40; magnoise=2.0/100.; phnoise=2*pi/180.; toll=le-10; tol2=le-10; 
% 

% Initialize 

DA=zeros (na , (nf+nr )) ; DAr=zeros (nx , (nf+nr) ) ; 
for n=l:nd 
k=fdof(n); 

DA(k,n)=l; 

k=frdof(n); 

DAr(k,n)=l; 

end 


% 

% Form Frequency-Response Vectors (Direct Solution) 
^AB,P]=Mdirect(W,MA,CA,KA,DA); 

% Simulate Noise 

jABn,Pn]=applynoise(AB,P,nsample,magnoise,phnoise); 

% Remove Response Vectors Corresponding to ”A” Structure 
AA= ABn (rdof , : ) ; 

% 

% Form Velocity and Displacement FRFs 
^.A,VA,XA,Hffj=calvxf(AA ) Pn ) Freq); Hrf=[|; 

% LHS of Equation 3.7 

VXFA=rVA;XA;-Hff ; 

% 

% Form Real/ Imaginary Partitions 
% RHS of Equation 3.7 
AA2=[real(AA) imag(AA)]; 

VXFA2= [real(VXFA) imag(VXFA)]; 

% 

% Step 1 - Identification of M"-1K and M'-lC and the resulting eigensolution 

% Solve Equation 3.7 

%CKDAE2=-AA2/VXFA2; 

gKDAE2,Vl ) Sl]=tls(VXFA2 , r AA2 , ,21,toll ) tol2,l); CKDAE2=CKDAE2’; 

% Complex Eigensolution Based on Estimated Minv*C and Minv*K 
CHA=CKDAE2(:,l:nx); 

KHA=CKDAE2(:,nx+l:ns); 

DHA=CKDAE2(: ,ns+ 1 :ns+nd) ; 

% Equation 3.11 
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ASAHATE=[CHA eye(nx);eye(nx) zeros (nx)]; 

BSAHATE= KHA zeros(nx); zeros(nx) -eye(nx)]; 

% Equation 3.12 Lambda=LambdaHat 
CTHETAE,CLAMAE]=eig(-BSAHATE,ASAHATE); 

CTHETAE,CLAMAE =eigensorti(CTHETAE,diag(CLAMAE)); 
zetass=-real (CL AM AE) . / abs(CLAMAE); 
omegss=imag(CLAMAE)./sqrt(l-zetass. " 2); 

% 

% Step 2-Identification of M, C, and K 

% 

% Estimation of Inverse Generalized Modal Parameters 

% 

CTX=CTHETAE( 1 :nx , : ) ; 

% 

E=zeros(nf*ne*nx,ns); R=zeros(ne,ns); Y=zeros(nf*ne*nx,l); 
for k=l:ne 
for r=l:ns 
% Equation 3.19 

R(k,r)=-W2(k)/(j*W(k)-CLAMAE(r)); 

end 

end 

% 

offset=ne*nx; 
for n=l:nf 
nof=n*offset-offset ; 
for k=l:ne 
kk=nf*k-(nf-n); 

DF=DA (: ,n) + DA(: ,nf+ 1 :nd) *Hrf (: ,kk) ; 
for r=l:ns 
% Equation 3.24 

E((l+(k-l)*nx)+nof:(k*nx)+nof,r)=R(k,r)*CTX(:,r)*CTX(:,r).’*DF; 

end 

for r=l:nx 

% LHS of Equation 3.23 
Y (nx* (k- 1 ) +r +nof ) = A A (r ,kk) ; 
end 
end 
end 
% 

% Inverse Generalized Mass Estimation 
% Solving Equation 3.23 
GASAIE=E\Y; 

% 

% Estimated Generalized Modal Parameter Vector - a r 
GASAE=1./GASAIE; 

% 



% Estimated Generalized Modal Parameter Vector - b r 
% Solving Equation 3.25 
GBSAE=-CLAMAE.*GASAE; 

% 

% Forming (Thetahat) *-l 
CTHI=inv(CTHETAE) ; 

% 

% Estimated State-variable A Matrix 
% Equation 3.26a 

ASAE=CTHI.’*diag(GASAE)*CTHI; 

% 

% Estimated State-variable B Matrix 
% Equation 3.26b 

BSAE=CTHI.’*diag(GBSAE)*CTHI; 

% 

% Extract System Matrices 
Ml =real(ASAE(l :nx,nx-|- 1 :ns)); 
M2=real(ASAE(nx+l:ns,l:nx)); 

M3=-real (BSAE(nx+ l:ns,nx-)-l :ns)) ; 

M=M3; 

C=real(ASAE(l:nx,l:nx)); 

K=real (BS AE( 1 :nx, 1 :nx)) ; 

G2=damping(M , K ,zetass( 1 :nx)) ; 
phie,eige]=eig(K,M); 

phie,eige]=eigensortr(phie,diag(eige)); freqe=sqrt(eige)/(2*pi); 


function [Haf ,H vf , Hxf ,Hff ] =cal vxf ( A A ,P, Freq) 

% [Haf,Hvf,Hff]=calvxf(AA,P,IYeq) 

% AA = Accelerations of the ”A” structure 
% P = applied forces 

% Haf = Aocelerance Frequency Response Function 


% Hvf = 
% Hxf = 
% Hff 


Hafl/Qw) 
Haf]/(-w'2) 


ne=lengthfFreq); nf=size(P,l); W=Freq*2*pi; W2=W. " 2; 
Haf=zerosfsizefAAll; Hvf=zeros(size(AA)); 

Hxf=zeros (size ( A A) ) ; Hff =zeros(size (P) ) ; 

% 

for k=l:ne 


nk2=nf*k; 

nkl=nk2-(nf-l); 

% 

% Form Autospectrum 
Pstar=conj(P(:,nkl:nk2)); 
Gff=P (: ,nkl :nk2) * Pstar; 
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Gffi=inv(Gff); 

Hff(:,nkl:nk2)=Gff*Gffi; 

% 

% Form [Hail/ (jw)-Mobility 
WB=eye(n/)*(l./G*W(k))); 

Hvf(:,nkl:nk2)=(AA(:,nkl:nk2)*WB)*Pstar*Gffi; 

% 

% Form [Hafl/(-w ~ 2)-Receptance 
W2B=-eye(nf)*l./W2(k); 

Hxf(: ,nkl :nk2)= ( A A (: ,nkl :nk2) * W2B) * Pstar*Gffi; 

% 

% Form Haf-Accelerance 
Haf(:,nkl:nk2)=AA(:,nkl:nk2)*Pstar*Gffi; 
end 



147 


function JABn,PnJ=applynoise(AB,P,nsample,magnoise,phnoise) 

% [ABn,Pn|=applynoise(AB,P,nsample,magnoise,phnoise) 

nx=size(AB,l); 

nf=size(P,l); 

ne=size(AB.2) /nf; 

ABn=zeros (size (A B) ) ; 

Pn=zeros(size(P)) ; 
for k=l:nx 
for kk=l:nf 

ABn(k,kk:nf:nf*ne)=addnoise(AB(k,kk:nf:nf*ne) ,nsample,magnoise,phnoise); 
end 
end 

for k=l:nf 

Pn(k,k:nf:nf*ne)=addnoise(P(k,k:nf:nf*ne),nsample,magnoise,phnoise); 

end 


function XN =addnoise(X , N avg,magnoise .phnoise) 
% XN=addnoise(X,Navg,magnoise,phnoise) 

% XN = output signal with the noise added 
% X = the original signal 
% Navg = the number of avergages 
% magnoise = the noise on the magnitude 
% phnoise = the noise on the phase 
[nx nf]=size(X); 

Xreal=real(X); 

Ximag=imag(X) ; 

% Compute magnitude and phase of base signal 
Xmag=sqrt( (Xreal). , '2 + (Ximag)."2 ); 
Xphase=atan2 (Ximag ,Xreal) ; 

% RMS value of base signal 
Hrms=norm(Xmag) /length(Xmag) ; 

% Percentage of RMS value of base signal 
noisescale=Hrms*magnoise; 

XN =zeros (size (X) 1 ; 

I=ones (size (Xmag) ) ; 
for k=l:Navg 
Nmag=zeros (size(Xmag)) ; 
Nphase=zeros(size(Xphase) ) ; 

% Calculate and add noise for magnitude 
Rl=rand(nx,nf) ; 


R2=rand(nx,nf); 
noiseM= (1-2 * R1 ) ; 
noiseM=noiseM-mean(noiseM) ; 
Nrms=norm(noiseM) /length (noiseM) ; 
SNmag=noisescale/Nrms; 
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N mag= S N mag* noiseM ; 

% Calculate and noise for phase 
noisePh= (1-2* R2) ; 
noisePh=noisePh-mean (noiseP h) ; 
Nphase=phnoise*noisePh; 

% 

XNmag=Xmag+Nmag; 

XNph=Xphase+Nphase; 

% 

XNreal=XNmag.*cos(XNph): 

XNimag=XNmag.*sin(XNph); 

XNnk=XNreal+j*XNimag; 

XN=XN+XNnk; 

end 

XN=XN/Navg; 



149 


function [X,V,S]=tls(A,B, rank, toll, tol2, comp) 

% Solves the Total Least Squares problem AX=B 
% pC, V,SJ=tls(A,B, rank, toll ,tol2, comp) 

% A Tne data matrix (mxn) 

% B The observation matrix (nxd) 

% rank The rank of [A;B1, if set to zero it will be computed 
% toll Two possibilities depending on the value of comp 
% Specifies the tolerance to be used when checking for multiplicity of 
% singular vales OR the estimated standard deviation of the error on 
% each element of the matrix C 

% tol2 Tolerance used to check for the singularity of the upper triangular 
% matrix F, see description 

% comp Specifies whether rank and/or toll is to be computed by the routine 
% comp = 1 rank is to be computed 

% =2 toll is to be computed (stdv is input by the user) 

% =3 neither toll nor rank is computed 

% =4 both rank and toll are to be computed (stdv is input) 

% X The solution to the TLS problem, the leading nxd part of the array 
% V The right singular vectors of [A;B] 

% S The singular values of C in descending order 


% Reference S.Van Huffel (ESAT Laboratory, KU Leuven). 
crank=l; ctol=l; 
if comp==3 
crank=0; 
ctol=0; 
else 

if comp==l 
ctol=0; 


end 

if comp==2 
crank=0; 

end 

end 

toll=max(Jtoll eps]); tol2=max([tol2 eps]); 
[m,n]=size(A); l=size(B,2); nl=n+l; 
A=zeros(n,l); k=max([m nil): p=min([m nl): 
fU,S,V]=sva([A,B],0); fc=diag(S); 


70 

% Step 2: Compute the rank approximation of [A B] 
smax=toll; 

% Compute smax if requested 
if ctol==l 


smax=sqrt (2. *k) *smax; 
end 

smax2=smax A 2; 
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% Compute the rank if requested 
if crank==l 
rank=p; 

while (rank>0) & (S(rank)<=smax) 
if rank > 0 
if S(rank) <= smax 
rank=rank-l; 
end 
end 
end 
end 

% Adjust the rank if S(rank) has multiplicity > 1 
flag=0; 

while flag==0 
rl=rank+l; 

while (rank>0) & ((S(rank) A 2-S(rl) "2)<=smax2) 
if rank > 0 

if (S(rank) * 2-S(rl ) A 2) <=smax2 
rank=rank-l; 

’Rank lowered because singular value had multiplicity > 1’ 
end 
end 
end 

rl=rank+l; 

% 

% Compute the Householder matrix Q and matrices F and Y 

% 

nll=max([n rl])+l; 

zero=0; 

i=nl; 

while ( ("zero) & (i>=nll) ) 
if ( ("zero) & (i>=nll) ) 
k=i-rank; 

WRK= V (i ,r 1 :rl +k-l) ; 

[WRK ,temp ,zero] =housh ( WRK ,k,tol2) ; 
if ("zero) 

V=tr2(V, WRK, temp,l ,i, rank, k); 
end 
i=i-l; 
end 
end 

% 

nl=n+l; 

if ((zero) | (abs(V(nl,nl))<=tol2)) 
rank=rank-l; 

‘Rank lowered because singular upper triangular matrix F’ 
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flag=0; 

else 

flag=l; 

end 

end 

X(l:n,l)=X(l:n,l)+(-l.)/V(nl,nl)*V(l:n,nl); 
for j=2:l 
nj=n+j; 
temp=V(nj,nj); 

ji=H; 

for i=l:n 

x (ij)=-(V(i,nj)+dot(V(nl:nl+jl-l,nj),X(i,l:jl)))/temp; 

end 

end 
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function [u,s,zero]=housh(uj,heps) 

%HOUSH Construct a householder transformation H=I-s*UU\ 

% 

% [U,S,ZERO] = HOUSH(U,J,Heps) 

% Constructs a Householder transformation H=I-s*UU’ that ‘mirrors’ a 
% vector u to the Jth unit vector. If NORM(U)<Eps then Zero=l [True] 

% 

% Reference: Adapted from “Computation of Zeros of Linear Multivariable 
% Systems”, A. Emami-Naeini, and P. Van Dooren; Automatica 

% Vol. 18, No. 4, pp. 415-430, 1982. 

s = sum(u.*u); 
alfa = sqrt(s); 

if (alfa<=heps), zero=l; return, end 

zero=0; 

dum = u(j); 

if dum>0, alfa=-aJfa; end 
uG) = uG)-alfa; 
s = 1 ./(s-alfa*dum); 

% Make u a column vector, 
u = u(:); 


function A=tr2(A,U,s,il,i2jl j2) 
inprod=0.0; 
for j=l:j2 

inprod=inprod-l-U G)*A(i jl-fj); 
end 

y=inprod*s; 
for j=l:j2 

A(ijl+j)=A(ijl+j)-UG)*y; 

end 

end 
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